xref: /dports/cad/gmsh/gmsh-4.9.2-source/api/gmsh.h_cwrap

// Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file in the Gmsh root directory for license information.
// Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.

#ifndef GMSH_H
#define GMSH_H

// This file redefines the Gmsh C++ API in terms of the C API (v4.9.2).
//
// This is provided as a convenience for users of the binary Gmsh SDK whose C++
// compiler ABI is not compatible with the ABI of the C++ compiler used to create
// the SDK (and who can thus not directly use the C++ API defined in `gmsh.h').
//
// To use this header file in your C++ code, simply rename it as `gmsh.h'.
//
// Note that using this header file will lead to (slightly) reduced performance
// compared to using the native Gmsh C++ API from the original `gmsh.h', as it
// entails additional data copies between this C++ wrapper, the C API and the
// native C++ code.
//
// Do not edit this file directly: it is automatically generated by `api/gen.py'.

#include <cmath>
#include <vector>
#include <string>
#include <utility>
#include <functional>

#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif

extern "C" {
  #include "gmshc.h"
}

namespace gmsh {

  // A geometrical entity in the Gmsh API is represented by two integers: its
  // dimension (dim = 0, 1, 2 or 3) and its tag (its unique, strictly positive
  // identifier). When dealing with multiple geometrical entities of possibly
  // different dimensions, the entities are packed as a vector of (dim, tag)
  // integer pairs.
  typedef std::vector<std::pair<int, int> > vectorpair;

}

namespace gmsh {

  template<typename t>
  inline void vector2ptr(const std::vector<t> &v, t **p, size_t *size)
  {
    *p = (t*)gmshMalloc(sizeof(t) * v.size());
    for(size_t i = 0; i < v.size(); ++i){
      (*p)[i] = v[i];
    }
    *size = v.size();
  }

  inline void vectorpair2intptr(const gmsh::vectorpair &v, int **p, size_t *size)
  {
    *p = (int*)gmshMalloc(sizeof(int) * v.size() * 2);
    for(size_t i = 0; i < v.size(); ++i){
      (*p)[i * 2 + 0] = v[i].first;
      (*p)[i * 2 + 1] = v[i].second;
    }
    *size = v.size() * 2;
  }

  inline void vectorstring2charptrptr(const std::vector<std::string> &v, char ***p, size_t *size)
  {
    *p = (char**)gmshMalloc(sizeof(char*) * v.size());
    for(size_t i = 0; i < v.size(); ++i){
      (*p)[i] = (char*)gmshMalloc(sizeof(char) * (v[i].size() + 1));
      for(size_t j = 0; j < v[i].size(); j++) (*p)[i][j] = v[i][j];
      (*p)[i][v[i].size()] = '\0';
    }
    *size = v.size();
  }

  template<typename t>
  inline void vectorvector2ptrptr(const std::vector<std::vector<t> > &v, t ***p, size_t **size, size_t *sizeSize)
  {
    *p = (t**)gmshMalloc(sizeof(t*) * v.size());
    *size = (size_t*)gmshMalloc(sizeof(size_t) * v.size());
    for(size_t i = 0; i < v.size(); ++i)
      vector2ptr(v[i], &((*p)[i]), &((*size)[i]));
    *sizeSize = v.size();
  }

  inline void throwLastError()
  {
     int ierr = 0;
     char *api_error_;
     gmshLoggerGetLastError(&api_error_, &ierr);
     if(ierr) throw "Could not get last error";
     std::string error = std::string(api_error_);
     gmshFree(api_error_);
     throw error;
  }

}

namespace gmsh { // Top-level functions

  // Initialize the Gmsh API. This must be called before any call to the other
  // functions in the API. If `argc' and `argv' (or just `argv' in Python or Julia)
  // are provided, they will be handled in the same way as the command line
  // arguments in the Gmsh app. If `readConfigFiles' is set, read system Gmsh
  // configuration files (gmshrc and gmsh-options). If `run' is set, run in the
  // same way as the Gmsh app, either interactively or in batch mode depending on
  // the command line arguments. If `run' is not set, initializing the API sets the
  // options "General.AbortOnError" to 2 and "General.Terminal" to 1. If compiled
  // with OpenMP support, it also sets the number of threads to
  // "General.NumThreads".
  inline void initialize(int argc = 0, char ** argv = 0,
                         const bool readConfigFiles = true,
                         const bool run = false)
  {
    int ierr = 0;
    gmshInitialize(argc, argv, (int)readConfigFiles, (int)run, &ierr);
    if(ierr) throwLastError();
  }

  // Finalize the Gmsh API. This must be called when you are done using the Gmsh
  // API.
  inline void finalize()
  {
    int ierr = 0;
    gmshFinalize(&ierr);
    if(ierr) throwLastError();
  }

  // Open a file. Equivalent to the `File->Open' menu in the Gmsh app. Handling of
  // the file depends on its extension and/or its contents: opening a file with
  // model data will create a new model.
  inline void open(const std::string & fileName)
  {
    int ierr = 0;
    gmshOpen(fileName.c_str(), &ierr);
    if(ierr) throwLastError();
  }

  // Merge a file. Equivalent to the `File->Merge' menu in the Gmsh app. Handling
  // of the file depends on its extension and/or its contents. Merging a file with
  // model data will add the data to the current model.
  inline void merge(const std::string & fileName)
  {
    int ierr = 0;
    gmshMerge(fileName.c_str(), &ierr);
    if(ierr) throwLastError();
  }

  // Write a file. The export format is determined by the file extension.
  inline void write(const std::string & fileName)
  {
    int ierr = 0;
    gmshWrite(fileName.c_str(), &ierr);
    if(ierr) throwLastError();
  }

  // Clear all loaded models and post-processing data, and add a new empty model.
  inline void clear()
  {
    int ierr = 0;
    gmshClear(&ierr);
    if(ierr) throwLastError();
  }

  namespace option { // Option handling functions

    // Set a numerical option to `value'. `name' is of the form "Category.Option"
    // or "Category[num].Option". Available categories and options are listed in
    // the Gmsh reference manual.
    inline void setNumber(const std::string & name,
                          const double value)
    {
      int ierr = 0;
      gmshOptionSetNumber(name.c_str(), value, &ierr);
      if(ierr) throwLastError();
    }

    // Get the `value' of a numerical option. `name' is of the form
    // "Category.Option" or "Category[num].Option". Available categories and
    // options are listed in the Gmsh reference manual.
    inline void getNumber(const std::string & name,
                          double & value)
    {
      int ierr = 0;
      gmshOptionGetNumber(name.c_str(), &value, &ierr);
      if(ierr) throwLastError();
    }

    // Set a string option to `value'. `name' is of the form "Category.Option" or
    // "Category[num].Option". Available categories and options are listed in the
    // Gmsh reference manual.
    inline void setString(const std::string & name,
                          const std::string & value)
    {
      int ierr = 0;
      gmshOptionSetString(name.c_str(), value.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Get the `value' of a string option. `name' is of the form "Category.Option"
    // or "Category[num].Option". Available categories and options are listed in
    // the Gmsh reference manual.
    inline void getString(const std::string & name,
                          std::string & value)
    {
      int ierr = 0;
      char *api_value_;
      gmshOptionGetString(name.c_str(), &api_value_, &ierr);
      if(ierr) throwLastError();
      value = std::string(api_value_); gmshFree(api_value_);
    }

    // Set a color option to the RGBA value (`r', `g', `b', `a'), where where `r',
    // `g', `b' and `a' should be integers between 0 and 255. `name' is of the form
    // "Category.Color.Option" or "Category[num].Color.Option". Available
    // categories and options are listed in the Gmsh reference manual. For
    // conciseness "Color." can be ommitted in `name'.
    inline void setColor(const std::string & name,
                         const int r,
                         const int g,
                         const int b,
                         const int a = 255)
    {
      int ierr = 0;
      gmshOptionSetColor(name.c_str(), r, g, b, a, &ierr);
      if(ierr) throwLastError();
    }

    // Get the `r', `g', `b', `a' value of a color option. `name' is of the form
    // "Category.Color.Option" or "Category[num].Color.Option". Available
    // categories and options are listed in the Gmsh reference manual. For
    // conciseness "Color." can be ommitted in `name'.
    inline void getColor(const std::string & name,
                         int & r,
                         int & g,
                         int & b,
                         int & a)
    {
      int ierr = 0;
      gmshOptionGetColor(name.c_str(), &r, &g, &b, &a, &ierr);
      if(ierr) throwLastError();
    }

  } // namespace option

  namespace model { // Model functions

    // Add a new model, with name `name', and set it as the current model.
    inline void add(const std::string & name)
    {
      int ierr = 0;
      gmshModelAdd(name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Remove the current model.
    inline void remove()
    {
      int ierr = 0;
      gmshModelRemove(&ierr);
      if(ierr) throwLastError();
    }

    // List the names of all models.
    inline void list(std::vector<std::string> & names)
    {
      int ierr = 0;
      char **api_names_; size_t api_names_n_;
      gmshModelList(&api_names_, &api_names_n_, &ierr);
      if(ierr) throwLastError();
      names.resize(api_names_n_); for(size_t i = 0; i < api_names_n_; ++i){ names[i] = std::string(api_names_[i]); gmshFree(api_names_[i]); } gmshFree(api_names_);
    }

    // Get the name of the current model.
    inline void getCurrent(std::string & name)
    {
      int ierr = 0;
      char *api_name_;
      gmshModelGetCurrent(&api_name_, &ierr);
      if(ierr) throwLastError();
      name = std::string(api_name_); gmshFree(api_name_);
    }

    // Set the current model to the model with name `name'. If several models have
    // the same name, select the one that was added first.
    inline void setCurrent(const std::string & name)
    {
      int ierr = 0;
      gmshModelSetCurrent(name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Get the file name (if any) associated with the current model. A file name is
    // associated when a model is read from a file on disk.
    inline void getFileName(std::string & fileName)
    {
      int ierr = 0;
      char *api_fileName_;
      gmshModelGetFileName(&api_fileName_, &ierr);
      if(ierr) throwLastError();
      fileName = std::string(api_fileName_); gmshFree(api_fileName_);
    }

    // Set the file name associated with the current model.
    inline void setFileName(const std::string & fileName)
    {
      int ierr = 0;
      gmshModelSetFileName(fileName.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Get all the entities in the current model. If `dim' is >= 0, return only the
    // entities of the specified dimension (e.g. points if `dim' == 0). The
    // entities are returned as a vector of (dim, tag) integer pairs.
    inline void getEntities(gmsh::vectorpair & dimTags,
                            const int dim = -1)
    {
      int ierr = 0;
      int *api_dimTags_; size_t api_dimTags_n_;
      gmshModelGetEntities(&api_dimTags_, &api_dimTags_n_, dim, &ierr);
      if(ierr) throwLastError();
      dimTags.resize(api_dimTags_n_ / 2); for(size_t i = 0; i < api_dimTags_n_ / 2; ++i){ dimTags[i].first = api_dimTags_[i * 2 + 0]; dimTags[i].second = api_dimTags_[i * 2 + 1]; } gmshFree(api_dimTags_);
    }

    // Set the name of the entity of dimension `dim' and tag `tag'.
    inline void setEntityName(const int dim,
                              const int tag,
                              const std::string & name)
    {
      int ierr = 0;
      gmshModelSetEntityName(dim, tag, name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Get the name of the entity of dimension `dim' and tag `tag'.
    inline void getEntityName(const int dim,
                              const int tag,
                              std::string & name)
    {
      int ierr = 0;
      char *api_name_;
      gmshModelGetEntityName(dim, tag, &api_name_, &ierr);
      if(ierr) throwLastError();
      name = std::string(api_name_); gmshFree(api_name_);
    }

    // Get all the physical groups in the current model. If `dim' is >= 0, return
    // only the entities of the specified dimension (e.g. physical points if `dim'
    // == 0). The entities are returned as a vector of (dim, tag) integer pairs.
    inline void getPhysicalGroups(gmsh::vectorpair & dimTags,
                                  const int dim = -1)
    {
      int ierr = 0;
      int *api_dimTags_; size_t api_dimTags_n_;
      gmshModelGetPhysicalGroups(&api_dimTags_, &api_dimTags_n_, dim, &ierr);
      if(ierr) throwLastError();
      dimTags.resize(api_dimTags_n_ / 2); for(size_t i = 0; i < api_dimTags_n_ / 2; ++i){ dimTags[i].first = api_dimTags_[i * 2 + 0]; dimTags[i].second = api_dimTags_[i * 2 + 1]; } gmshFree(api_dimTags_);
    }

    // Get the tags of the model entities making up the physical group of dimension
    // `dim' and tag `tag'.
    inline void getEntitiesForPhysicalGroup(const int dim,
                                            const int tag,
                                            std::vector<int> & tags)
    {
      int ierr = 0;
      int *api_tags_; size_t api_tags_n_;
      gmshModelGetEntitiesForPhysicalGroup(dim, tag, &api_tags_, &api_tags_n_, &ierr);
      if(ierr) throwLastError();
      tags.assign(api_tags_, api_tags_ + api_tags_n_); gmshFree(api_tags_);
    }

    // Get the tags of the physical groups (if any) to which the model entity of
    // dimension `dim' and tag `tag' belongs.
    inline void getPhysicalGroupsForEntity(const int dim,
                                           const int tag,
                                           std::vector<int> & physicalTags)
    {
      int ierr = 0;
      int *api_physicalTags_; size_t api_physicalTags_n_;
      gmshModelGetPhysicalGroupsForEntity(dim, tag, &api_physicalTags_, &api_physicalTags_n_, &ierr);
      if(ierr) throwLastError();
      physicalTags.assign(api_physicalTags_, api_physicalTags_ + api_physicalTags_n_); gmshFree(api_physicalTags_);
    }

    // Add a physical group of dimension `dim', grouping the model entities with
    // tags `tags'. Return the tag of the physical group, equal to `tag' if `tag'
    // is positive, or a new tag if `tag' < 0.
    inline int addPhysicalGroup(const int dim,
                                const std::vector<int> & tags,
                                const int tag = -1)
    {
      int ierr = 0;
      int *api_tags_; size_t api_tags_n_; vector2ptr(tags, &api_tags_, &api_tags_n_);
      int result_api_ = gmshModelAddPhysicalGroup(dim, api_tags_, api_tags_n_, tag, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_tags_);
      return result_api_;
    }

    // Remove the physical groups `dimTags' from the current model. If `dimTags' is
    // empty, remove all groups.
    inline void removePhysicalGroups(const gmsh::vectorpair & dimTags = gmsh::vectorpair())
    {
      int ierr = 0;
      int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
      gmshModelRemovePhysicalGroups(api_dimTags_, api_dimTags_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_dimTags_);
    }

    // Set the name of the physical group of dimension `dim' and tag `tag'.
    inline void setPhysicalName(const int dim,
                                const int tag,
                                const std::string & name)
    {
      int ierr = 0;
      gmshModelSetPhysicalName(dim, tag, name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Remove the physical name `name' from the current model.
    inline void removePhysicalName(const std::string & name)
    {
      int ierr = 0;
      gmshModelRemovePhysicalName(name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Get the name of the physical group of dimension `dim' and tag `tag'.
    inline void getPhysicalName(const int dim,
                                const int tag,
                                std::string & name)
    {
      int ierr = 0;
      char *api_name_;
      gmshModelGetPhysicalName(dim, tag, &api_name_, &ierr);
      if(ierr) throwLastError();
      name = std::string(api_name_); gmshFree(api_name_);
    }

    // Set the tag of the entity of dimension `dim' and tag `tag' to the new value
    // `newTag'.
    inline void setTag(const int dim,
                       const int tag,
                       const int newTag)
    {
      int ierr = 0;
      gmshModelSetTag(dim, tag, newTag, &ierr);
      if(ierr) throwLastError();
    }

    // Get the boundary of the model entities `dimTags'. Return in `outDimTags' the
    // boundary of the individual entities (if `combined' is false) or the boundary
    // of the combined geometrical shape formed by all input entities (if
    // `combined' is true). Return tags multiplied by the sign of the boundary
    // entity if `oriented' is true. Apply the boundary operator recursively down
    // to dimension 0 (i.e. to points) if `recursive' is true.
    inline void getBoundary(const gmsh::vectorpair & dimTags,
                            gmsh::vectorpair & outDimTags,
                            const bool combined = true,
                            const bool oriented = true,
                            const bool recursive = false)
    {
      int ierr = 0;
      int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
      int *api_outDimTags_; size_t api_outDimTags_n_;
      gmshModelGetBoundary(api_dimTags_, api_dimTags_n_, &api_outDimTags_, &api_outDimTags_n_, (int)combined, (int)oriented, (int)recursive, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_dimTags_);
      outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
    }

    // Get the upward and downward adjacencies of the model entity of dimension
    // `dim' and tag `tag'. The `upward' vector returns the adjacent entities of
    // dimension `dim' + 1; the `downward' vector returns the adjacent entities of
    // dimension `dim' - 1.
    inline void getAdjacencies(const int dim,
                               const int tag,
                               std::vector<int> & upward,
                               std::vector<int> & downward)
    {
      int ierr = 0;
      int *api_upward_; size_t api_upward_n_;
      int *api_downward_; size_t api_downward_n_;
      gmshModelGetAdjacencies(dim, tag, &api_upward_, &api_upward_n_, &api_downward_, &api_downward_n_, &ierr);
      if(ierr) throwLastError();
      upward.assign(api_upward_, api_upward_ + api_upward_n_); gmshFree(api_upward_);
      downward.assign(api_downward_, api_downward_ + api_downward_n_); gmshFree(api_downward_);
    }

    // Get the model entities in the bounding box defined by the two points
    // (`xmin', `ymin', `zmin') and (`xmax', `ymax', `zmax'). If `dim' is >= 0,
    // return only the entities of the specified dimension (e.g. points if `dim' ==
    // 0).
    inline void getEntitiesInBoundingBox(const double xmin,
                                         const double ymin,
                                         const double zmin,
                                         const double xmax,
                                         const double ymax,
                                         const double zmax,
                                         gmsh::vectorpair & tags,
                                         const int dim = -1)
    {
      int ierr = 0;
      int *api_tags_; size_t api_tags_n_;
      gmshModelGetEntitiesInBoundingBox(xmin, ymin, zmin, xmax, ymax, zmax, &api_tags_, &api_tags_n_, dim, &ierr);
      if(ierr) throwLastError();
      tags.resize(api_tags_n_ / 2); for(size_t i = 0; i < api_tags_n_ / 2; ++i){ tags[i].first = api_tags_[i * 2 + 0]; tags[i].second = api_tags_[i * 2 + 1]; } gmshFree(api_tags_);
    }

    // Get the bounding box (`xmin', `ymin', `zmin'), (`xmax', `ymax', `zmax') of
    // the model entity of dimension `dim' and tag `tag'. If `dim' and `tag' are
    // negative, get the bounding box of the whole model.
    inline void getBoundingBox(const int dim,
                               const int tag,
                               double & xmin,
                               double & ymin,
                               double & zmin,
                               double & xmax,
                               double & ymax,
                               double & zmax)
    {
      int ierr = 0;
      gmshModelGetBoundingBox(dim, tag, &xmin, &ymin, &zmin, &xmax, &ymax, &zmax, &ierr);
      if(ierr) throwLastError();
    }

    // Return the geometrical dimension of the current model.
    inline int getDimension()
    {
      int ierr = 0;
      int result_api_ = gmshModelGetDimension(&ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

    // Add a discrete model entity (defined by a mesh) of dimension `dim' in the
    // current model. Return the tag of the new discrete entity, equal to `tag' if
    // `tag' is positive, or a new tag if `tag' < 0. `boundary' specifies the tags
    // of the entities on the boundary of the discrete entity, if any. Specifying
    // `boundary' allows Gmsh to construct the topology of the overall model.
    inline int addDiscreteEntity(const int dim,
                                 const int tag = -1,
                                 const std::vector<int> & boundary = std::vector<int>())
    {
      int ierr = 0;
      int *api_boundary_; size_t api_boundary_n_; vector2ptr(boundary, &api_boundary_, &api_boundary_n_);
      int result_api_ = gmshModelAddDiscreteEntity(dim, tag, api_boundary_, api_boundary_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_boundary_);
      return result_api_;
    }

    // Remove the entities `dimTags' of the current model, provided that they are
    // not on the boundary of (or embedded in) higher-dimensional entities. If
    // `recursive' is true, remove all the entities on their boundaries, down to
    // dimension 0.
    inline void removeEntities(const gmsh::vectorpair & dimTags,
                               const bool recursive = false)
    {
      int ierr = 0;
      int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
      gmshModelRemoveEntities(api_dimTags_, api_dimTags_n_, (int)recursive, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_dimTags_);
    }

    // Remove the entity name `name' from the current model.
    inline void removeEntityName(const std::string & name)
    {
      int ierr = 0;
      gmshModelRemoveEntityName(name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Get the type of the entity of dimension `dim' and tag `tag'.
    inline void getType(const int dim,
                        const int tag,
                        std::string & entityType)
    {
      int ierr = 0;
      char *api_entityType_;
      gmshModelGetType(dim, tag, &api_entityType_, &ierr);
      if(ierr) throwLastError();
      entityType = std::string(api_entityType_); gmshFree(api_entityType_);
    }

    // In a partitioned model, get the parent of the entity of dimension `dim' and
    // tag `tag', i.e. from which the entity is a part of, if any. `parentDim' and
    // `parentTag' are set to -1 if the entity has no parent.
    inline void getParent(const int dim,
                          const int tag,
                          int & parentDim,
                          int & parentTag)
    {
      int ierr = 0;
      gmshModelGetParent(dim, tag, &parentDim, &parentTag, &ierr);
      if(ierr) throwLastError();
    }

    // Return the number of partitions in the model.
    inline int getNumberOfPartitions()
    {
      int ierr = 0;
      int result_api_ = gmshModelGetNumberOfPartitions(&ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

    // In a partitioned model, return the tags of the partition(s) to which the
    // entity belongs.
    inline void getPartitions(const int dim,
                              const int tag,
                              std::vector<int> & partitions)
    {
      int ierr = 0;
      int *api_partitions_; size_t api_partitions_n_;
      gmshModelGetPartitions(dim, tag, &api_partitions_, &api_partitions_n_, &ierr);
      if(ierr) throwLastError();
      partitions.assign(api_partitions_, api_partitions_ + api_partitions_n_); gmshFree(api_partitions_);
    }

    // Evaluate the parametrization of the entity of dimension `dim' and tag `tag'
    // at the parametric coordinates `parametricCoord'. Only valid for `dim' equal
    // to 0 (with empty `parametricCoord'), 1 (with `parametricCoord' containing
    // parametric coordinates on the curve) or 2 (with `parametricCoord' containing
    // pairs of u, v parametric coordinates on the surface, concatenated: [p1u,
    // p1v, p2u, ...]). Return triplets of x, y, z coordinates in `coord',
    // concatenated: [p1x, p1y, p1z, p2x, ...].
    inline void getValue(const int dim,
                         const int tag,
                         const std::vector<double> & parametricCoord,
                         std::vector<double> & coord)
    {
      int ierr = 0;
      double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
      double *api_coord_; size_t api_coord_n_;
      gmshModelGetValue(dim, tag, api_parametricCoord_, api_parametricCoord_n_, &api_coord_, &api_coord_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_parametricCoord_);
      coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
    }

    // Evaluate the derivative of the parametrization of the entity of dimension
    // `dim' and tag `tag' at the parametric coordinates `parametricCoord'. Only
    // valid for `dim' equal to 1 (with `parametricCoord' containing parametric
    // coordinates on the curve) or 2 (with `parametricCoord' containing pairs of
    // u, v parametric coordinates on the surface, concatenated: [p1u, p1v, p2u,
    // ...]). For `dim' equal to 1 return the x, y, z components of the derivative
    // with respect to u [d1ux, d1uy, d1uz, d2ux, ...]; for `dim' equal to 2 return
    // the x, y, z components of the derivative with respect to u and v: [d1ux,
    // d1uy, d1uz, d1vx, d1vy, d1vz, d2ux, ...].
    inline void getDerivative(const int dim,
                              const int tag,
                              const std::vector<double> & parametricCoord,
                              std::vector<double> & derivatives)
    {
      int ierr = 0;
      double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
      double *api_derivatives_; size_t api_derivatives_n_;
      gmshModelGetDerivative(dim, tag, api_parametricCoord_, api_parametricCoord_n_, &api_derivatives_, &api_derivatives_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_parametricCoord_);
      derivatives.assign(api_derivatives_, api_derivatives_ + api_derivatives_n_); gmshFree(api_derivatives_);
    }

    // Evaluate the second derivative of the parametrization of the entity of
    // dimension `dim' and tag `tag' at the parametric coordinates
    // `parametricCoord'. Only valid for `dim' equal to 1 (with `parametricCoord'
    // containing parametric coordinates on the curve) or 2 (with `parametricCoord'
    // containing pairs of u, v parametric coordinates on the surface,
    // concatenated: [p1u, p1v, p2u, ...]). For `dim' equal to 1 return the x, y, z
    // components of the second derivative with respect to u [d1uux, d1uuy, d1uuz,
    // d2uux, ...]; for `dim' equal to 2 return the x, y, z components of the
    // second derivative with respect to u and v, and the mixed derivative with
    // respect to u and v: [d1uux, d1uuy, d1uuz, d1vvx, d1vvy, d1vvz, d1uvx, d1uvy,
    // d1uvz, d2uux, ...].
    inline void getSecondDerivative(const int dim,
                                    const int tag,
                                    const std::vector<double> & parametricCoord,
                                    std::vector<double> & derivatives)
    {
      int ierr = 0;
      double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
      double *api_derivatives_; size_t api_derivatives_n_;
      gmshModelGetSecondDerivative(dim, tag, api_parametricCoord_, api_parametricCoord_n_, &api_derivatives_, &api_derivatives_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_parametricCoord_);
      derivatives.assign(api_derivatives_, api_derivatives_ + api_derivatives_n_); gmshFree(api_derivatives_);
    }

    // Evaluate the (maximum) curvature of the entity of dimension `dim' and tag
    // `tag' at the parametric coordinates `parametricCoord'. Only valid for `dim'
    // equal to 1 (with `parametricCoord' containing parametric coordinates on the
    // curve) or 2 (with `parametricCoord' containing pairs of u, v parametric
    // coordinates on the surface, concatenated: [p1u, p1v, p2u, ...]).
    inline void getCurvature(const int dim,
                             const int tag,
                             const std::vector<double> & parametricCoord,
                             std::vector<double> & curvatures)
    {
      int ierr = 0;
      double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
      double *api_curvatures_; size_t api_curvatures_n_;
      gmshModelGetCurvature(dim, tag, api_parametricCoord_, api_parametricCoord_n_, &api_curvatures_, &api_curvatures_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_parametricCoord_);
      curvatures.assign(api_curvatures_, api_curvatures_ + api_curvatures_n_); gmshFree(api_curvatures_);
    }

    // Evaluate the principal curvatures of the surface with tag `tag' at the
    // parametric coordinates `parametricCoord', as well as their respective
    // directions. `parametricCoord' are given by pair of u and v coordinates,
    // concatenated: [p1u, p1v, p2u, ...].
    inline void getPrincipalCurvatures(const int tag,
                                       const std::vector<double> & parametricCoord,
                                       std::vector<double> & curvatureMax,
                                       std::vector<double> & curvatureMin,
                                       std::vector<double> & directionMax,
                                       std::vector<double> & directionMin)
    {
      int ierr = 0;
      double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
      double *api_curvatureMax_; size_t api_curvatureMax_n_;
      double *api_curvatureMin_; size_t api_curvatureMin_n_;
      double *api_directionMax_; size_t api_directionMax_n_;
      double *api_directionMin_; size_t api_directionMin_n_;
      gmshModelGetPrincipalCurvatures(tag, api_parametricCoord_, api_parametricCoord_n_, &api_curvatureMax_, &api_curvatureMax_n_, &api_curvatureMin_, &api_curvatureMin_n_, &api_directionMax_, &api_directionMax_n_, &api_directionMin_, &api_directionMin_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_parametricCoord_);
      curvatureMax.assign(api_curvatureMax_, api_curvatureMax_ + api_curvatureMax_n_); gmshFree(api_curvatureMax_);
      curvatureMin.assign(api_curvatureMin_, api_curvatureMin_ + api_curvatureMin_n_); gmshFree(api_curvatureMin_);
      directionMax.assign(api_directionMax_, api_directionMax_ + api_directionMax_n_); gmshFree(api_directionMax_);
      directionMin.assign(api_directionMin_, api_directionMin_ + api_directionMin_n_); gmshFree(api_directionMin_);
    }

    // Get the normal to the surface with tag `tag' at the parametric coordinates
    // `parametricCoord'. `parametricCoord' are given by pairs of u and v
    // coordinates, concatenated: [p1u, p1v, p2u, ...]. `normals' are returned as
    // triplets of x, y, z components, concatenated: [n1x, n1y, n1z, n2x, ...].
    inline void getNormal(const int tag,
                          const std::vector<double> & parametricCoord,
                          std::vector<double> & normals)
    {
      int ierr = 0;
      double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
      double *api_normals_; size_t api_normals_n_;
      gmshModelGetNormal(tag, api_parametricCoord_, api_parametricCoord_n_, &api_normals_, &api_normals_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_parametricCoord_);
      normals.assign(api_normals_, api_normals_ + api_normals_n_); gmshFree(api_normals_);
    }

    // Get the parametric coordinates `parametricCoord' for the points `coord' on
    // the entity of dimension `dim' and tag `tag'. `coord' are given as triplets
    // of x, y, z coordinates, concatenated: [p1x, p1y, p1z, p2x, ...].
    // `parametricCoord' returns the parametric coordinates t on the curve (if
    // `dim' = 1) or pairs of u and v coordinates concatenated on the surface (if
    // `dim' = 2), i.e. [p1t, p2t, ...] or [p1u, p1v, p2u, ...].
    inline void getParametrization(const int dim,
                                   const int tag,
                                   const std::vector<double> & coord,
                                   std::vector<double> & parametricCoord)
    {
      int ierr = 0;
      double *api_coord_; size_t api_coord_n_; vector2ptr(coord, &api_coord_, &api_coord_n_);
      double *api_parametricCoord_; size_t api_parametricCoord_n_;
      gmshModelGetParametrization(dim, tag, api_coord_, api_coord_n_, &api_parametricCoord_, &api_parametricCoord_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_coord_);
      parametricCoord.assign(api_parametricCoord_, api_parametricCoord_ + api_parametricCoord_n_); gmshFree(api_parametricCoord_);
    }

    // Get the `min' and `max' bounds of the parametric coordinates for the entity
    // of dimension `dim' and tag `tag'.
    inline void getParametrizationBounds(const int dim,
                                         const int tag,
                                         std::vector<double> & min,
                                         std::vector<double> & max)
    {
      int ierr = 0;
      double *api_min_; size_t api_min_n_;
      double *api_max_; size_t api_max_n_;
      gmshModelGetParametrizationBounds(dim, tag, &api_min_, &api_min_n_, &api_max_, &api_max_n_, &ierr);
      if(ierr) throwLastError();
      min.assign(api_min_, api_min_ + api_min_n_); gmshFree(api_min_);
      max.assign(api_max_, api_max_ + api_max_n_); gmshFree(api_max_);
    }

    // Check if the coordinates (or the parametric coordinates if `parametric' is
    // set) provided in `coord' correspond to points inside the entity of dimension
    // `dim' and tag `tag', and return the number of points inside. This feature is
    // only available for a subset of entities, depending on the underyling
    // geometrical representation.
    inline int isInside(const int dim,
                        const int tag,
                        const std::vector<double> & coord,
                        const bool parametric = false)
    {
      int ierr = 0;
      double *api_coord_; size_t api_coord_n_; vector2ptr(coord, &api_coord_, &api_coord_n_);
      int result_api_ = gmshModelIsInside(dim, tag, api_coord_, api_coord_n_, (int)parametric, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_coord_);
      return result_api_;
    }

    // Get the points `closestCoord' on the entity of dimension `dim' and tag `tag'
    // to the points `coord', by orthogonal projection. `coord' and `closestCoord'
    // are given as triplets of x, y, z coordinates, concatenated: [p1x, p1y, p1z,
    // p2x, ...]. `parametricCoord' returns the parametric coordinates t on the
    // curve (if `dim' = 1) or pairs of u and v coordinates concatenated on the
    // surface (if `dim' = 2), i.e. [p1t, p2t, ...] or [p1u, p1v, p2u, ...].
    inline void getClosestPoint(const int dim,
                                const int tag,
                                const std::vector<double> & coord,
                                std::vector<double> & closestCoord,
                                std::vector<double> & parametricCoord)
    {
      int ierr = 0;
      double *api_coord_; size_t api_coord_n_; vector2ptr(coord, &api_coord_, &api_coord_n_);
      double *api_closestCoord_; size_t api_closestCoord_n_;
      double *api_parametricCoord_; size_t api_parametricCoord_n_;
      gmshModelGetClosestPoint(dim, tag, api_coord_, api_coord_n_, &api_closestCoord_, &api_closestCoord_n_, &api_parametricCoord_, &api_parametricCoord_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_coord_);
      closestCoord.assign(api_closestCoord_, api_closestCoord_ + api_closestCoord_n_); gmshFree(api_closestCoord_);
      parametricCoord.assign(api_parametricCoord_, api_parametricCoord_ + api_parametricCoord_n_); gmshFree(api_parametricCoord_);
    }

    // Reparametrize the boundary entity (point or curve, i.e. with `dim' == 0 or
    // `dim' == 1) of tag `tag' on the surface `surfaceTag'. If `dim' == 1,
    // reparametrize all the points corresponding to the parametric coordinates
    // `parametricCoord'. Multiple matches in case of periodic surfaces can be
    // selected with `which'. This feature is only available for a subset of
    // entities, depending on the underyling geometrical representation.
    inline void reparametrizeOnSurface(const int dim,
                                       const int tag,
                                       const std::vector<double> & parametricCoord,
                                       const int surfaceTag,
                                       std::vector<double> & surfaceParametricCoord,
                                       const int which = 0)
    {
      int ierr = 0;
      double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
      double *api_surfaceParametricCoord_; size_t api_surfaceParametricCoord_n_;
      gmshModelReparametrizeOnSurface(dim, tag, api_parametricCoord_, api_parametricCoord_n_, surfaceTag, &api_surfaceParametricCoord_, &api_surfaceParametricCoord_n_, which, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_parametricCoord_);
      surfaceParametricCoord.assign(api_surfaceParametricCoord_, api_surfaceParametricCoord_ + api_surfaceParametricCoord_n_); gmshFree(api_surfaceParametricCoord_);
    }

    // Set the visibility of the model entities `dimTags' to `value'. Apply the
    // visibility setting recursively if `recursive' is true.
    inline void setVisibility(const gmsh::vectorpair & dimTags,
                              const int value,
                              const bool recursive = false)
    {
      int ierr = 0;
      int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
      gmshModelSetVisibility(api_dimTags_, api_dimTags_n_, value, (int)recursive, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_dimTags_);
    }

    // Get the visibility of the model entity of dimension `dim' and tag `tag'.
    inline void getVisibility(const int dim,
                              const int tag,
                              int & value)
    {
      int ierr = 0;
      gmshModelGetVisibility(dim, tag, &value, &ierr);
      if(ierr) throwLastError();
    }

    // Set the global visibility of the model per window to `value', where
    // `windowIndex' identifies the window in the window list.
    inline void setVisibilityPerWindow(const int value,
                                       const int windowIndex = 0)
    {
      int ierr = 0;
      gmshModelSetVisibilityPerWindow(value, windowIndex, &ierr);
      if(ierr) throwLastError();
    }

    // Set the color of the model entities `dimTags' to the RGBA value (`r', `g',
    // `b', `a'), where `r', `g', `b' and `a' should be integers between 0 and 255.
    // Apply the color setting recursively if `recursive' is true.
    inline void setColor(const gmsh::vectorpair & dimTags,
                         const int r,
                         const int g,
                         const int b,
                         const int a = 255,
                         const bool recursive = false)
    {
      int ierr = 0;
      int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
      gmshModelSetColor(api_dimTags_, api_dimTags_n_, r, g, b, a, (int)recursive, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_dimTags_);
    }

    // Get the color of the model entity of dimension `dim' and tag `tag'.
    inline void getColor(const int dim,
                         const int tag,
                         int & r,
                         int & g,
                         int & b,
                         int & a)
    {
      int ierr = 0;
      gmshModelGetColor(dim, tag, &r, &g, &b, &a, &ierr);
      if(ierr) throwLastError();
    }

    // Set the `x', `y', `z' coordinates of a geometrical point.
    inline void setCoordinates(const int tag,
                               const double x,
                               const double y,
                               const double z)
    {
      int ierr = 0;
      gmshModelSetCoordinates(tag, x, y, z, &ierr);
      if(ierr) throwLastError();
    }

    namespace mesh { // Mesh functions

      // Generate a mesh of the current model, up to dimension `dim' (0, 1, 2 or
      // 3).
      inline void generate(const int dim = 3)
      {
        int ierr = 0;
        gmshModelMeshGenerate(dim, &ierr);
        if(ierr) throwLastError();
      }

      // Partition the mesh of the current model into `numPart' partitions.
      // Optionally, `elementTags' and `partitions' can be provided to specify the
      // partition of each element explicitly.
      inline void partition(const int numPart,
                            const std::vector<std::size_t> & elementTags = std::vector<std::size_t>(),
                            const std::vector<int> & partitions = std::vector<int>())
      {
        int ierr = 0;
        size_t *api_elementTags_; size_t api_elementTags_n_; vector2ptr(elementTags, &api_elementTags_, &api_elementTags_n_);
        int *api_partitions_; size_t api_partitions_n_; vector2ptr(partitions, &api_partitions_, &api_partitions_n_);
        gmshModelMeshPartition(numPart, api_elementTags_, api_elementTags_n_, api_partitions_, api_partitions_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_elementTags_);
        gmshFree(api_partitions_);
      }

      // Unpartition the mesh of the current model.
      inline void unpartition()
      {
        int ierr = 0;
        gmshModelMeshUnpartition(&ierr);
        if(ierr) throwLastError();
      }

      // Optimize the mesh of the current model using `method' (empty for default
      // tetrahedral mesh optimizer, "Netgen" for Netgen optimizer, "HighOrder" for
      // direct high-order mesh optimizer, "HighOrderElastic" for high-order
      // elastic smoother, "HighOrderFastCurving" for fast curving algorithm,
      // "Laplace2D" for Laplace smoothing, "Relocate2D" and "Relocate3D" for node
      // relocation). If `force' is set apply the optimization also to discrete
      // entities. If `dimTags' is given, only apply the optimizer to the given
      // entities.
      inline void optimize(const std::string & method,
                           const bool force = false,
                           const int niter = 1,
                           const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshOptimize(method.c_str(), (int)force, niter, api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Recombine the mesh of the current model.
      inline void recombine()
      {
        int ierr = 0;
        gmshModelMeshRecombine(&ierr);
        if(ierr) throwLastError();
      }

      // Refine the mesh of the current model by uniformly splitting the elements.
      inline void refine()
      {
        int ierr = 0;
        gmshModelMeshRefine(&ierr);
        if(ierr) throwLastError();
      }

      // Set the order of the elements in the mesh of the current model to `order'.
      inline void setOrder(const int order)
      {
        int ierr = 0;
        gmshModelMeshSetOrder(order, &ierr);
        if(ierr) throwLastError();
      }

      // Get the last entities (if any) where a meshing error occurred. Currently
      // only populated by the new 3D meshing algorithms.
      inline void getLastEntityError(gmsh::vectorpair & dimTags)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_;
        gmshModelMeshGetLastEntityError(&api_dimTags_, &api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        dimTags.resize(api_dimTags_n_ / 2); for(size_t i = 0; i < api_dimTags_n_ / 2; ++i){ dimTags[i].first = api_dimTags_[i * 2 + 0]; dimTags[i].second = api_dimTags_[i * 2 + 1]; } gmshFree(api_dimTags_);
      }

      // Get the last nodes (if any) where a meshing error occurred. Currently only
      // populated by the new 3D meshing algorithms.
      inline void getLastNodeError(std::vector<std::size_t> & nodeTags)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        gmshModelMeshGetLastNodeError(&api_nodeTags_, &api_nodeTags_n_, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
      }

      // Clear the mesh, i.e. delete all the nodes and elements, for the entities
      // `dimTags'. If `dimTags' is empty, clear the whole mesh. Note that the mesh
      // of an entity can only be cleared if this entity is not on the boundary of
      // another entity with a non-empty mesh.
      inline void clear(const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshClear(api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Reverse the orientation of the elements in the entities `dimTags'. If
      // `dimTags' is empty, reverse the orientation of the elements in the whole
      // mesh.
      inline void reverse(const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshReverse(api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Apply the affine transformation `affineTransform' (16 entries of a 4x4
      // matrix, by row; only the 12 first can be provided for convenience) to the
      // coordinates of the nodes classified on the entities `dimTags'. If
      // `dimTags' is empty, transform all the nodes in the mesh.
      inline void affineTransform(const std::vector<double> & affineTransform,
                                  const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        double *api_affineTransform_; size_t api_affineTransform_n_; vector2ptr(affineTransform, &api_affineTransform_, &api_affineTransform_n_);
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshAffineTransform(api_affineTransform_, api_affineTransform_n_, api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_affineTransform_);
        gmshFree(api_dimTags_);
      }

      // Get the nodes classified on the entity of dimension `dim' and tag `tag'.
      // If `tag' < 0, get the nodes for all entities of dimension `dim'. If `dim'
      // and `tag' are negative, get all the nodes in the mesh. `nodeTags' contains
      // the node tags (their unique, strictly positive identification numbers).
      // `coord' is a vector of length 3 times the length of `nodeTags' that
      // contains the x, y, z coordinates of the nodes, concatenated: [n1x, n1y,
      // n1z, n2x, ...]. If `dim' >= 0 and `returnParamtricCoord' is set,
      // `parametricCoord' contains the parametric coordinates ([u1, u2, ...] or
      // [u1, v1, u2, ...]) of the nodes, if available. The length of
      // `parametricCoord' can be 0 or `dim' times the length of `nodeTags'. If
      // `includeBoundary' is set, also return the nodes classified on the boundary
      // of the entity (which will be reparametrized on the entity if `dim' >= 0 in
      // order to compute their parametric coordinates).
      inline void getNodes(std::vector<std::size_t> & nodeTags,
                           std::vector<double> & coord,
                           std::vector<double> & parametricCoord,
                           const int dim = -1,
                           const int tag = -1,
                           const bool includeBoundary = false,
                           const bool returnParametricCoord = true)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        double *api_coord_; size_t api_coord_n_;
        double *api_parametricCoord_; size_t api_parametricCoord_n_;
        gmshModelMeshGetNodes(&api_nodeTags_, &api_nodeTags_n_, &api_coord_, &api_coord_n_, &api_parametricCoord_, &api_parametricCoord_n_, dim, tag, (int)includeBoundary, (int)returnParametricCoord, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
        parametricCoord.assign(api_parametricCoord_, api_parametricCoord_ + api_parametricCoord_n_); gmshFree(api_parametricCoord_);
      }

      // Get the nodes classified on the entity of tag `tag', for all the elements
      // of type `elementType'. The other arguments are treated as in `getNodes'.
      inline void getNodesByElementType(const int elementType,
                                        std::vector<std::size_t> & nodeTags,
                                        std::vector<double> & coord,
                                        std::vector<double> & parametricCoord,
                                        const int tag = -1,
                                        const bool returnParametricCoord = true)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        double *api_coord_; size_t api_coord_n_;
        double *api_parametricCoord_; size_t api_parametricCoord_n_;
        gmshModelMeshGetNodesByElementType(elementType, &api_nodeTags_, &api_nodeTags_n_, &api_coord_, &api_coord_n_, &api_parametricCoord_, &api_parametricCoord_n_, tag, (int)returnParametricCoord, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
        parametricCoord.assign(api_parametricCoord_, api_parametricCoord_ + api_parametricCoord_n_); gmshFree(api_parametricCoord_);
      }

      // Get the coordinates and the parametric coordinates (if any) of the node
      // with tag `tag', as well as the dimension `dim' and tag `tag' of the entity
      // on which the node is classified. This function relies on an internal cache
      // (a vector in case of dense node numbering, a map otherwise); for large
      // meshes accessing nodes in bulk is often preferable.
      inline void getNode(const std::size_t nodeTag,
                          std::vector<double> & coord,
                          std::vector<double> & parametricCoord,
                          int & dim,
                          int & tag)
      {
        int ierr = 0;
        double *api_coord_; size_t api_coord_n_;
        double *api_parametricCoord_; size_t api_parametricCoord_n_;
        gmshModelMeshGetNode(nodeTag, &api_coord_, &api_coord_n_, &api_parametricCoord_, &api_parametricCoord_n_, &dim, &tag, &ierr);
        if(ierr) throwLastError();
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
        parametricCoord.assign(api_parametricCoord_, api_parametricCoord_ + api_parametricCoord_n_); gmshFree(api_parametricCoord_);
      }

      // Set the coordinates and the parametric coordinates (if any) of the node
      // with tag `tag'. This function relies on an internal cache (a vector in
      // case of dense node numbering, a map otherwise); for large meshes accessing
      // nodes in bulk is often preferable.
      inline void setNode(const std::size_t nodeTag,
                          const std::vector<double> & coord,
                          const std::vector<double> & parametricCoord)
      {
        int ierr = 0;
        double *api_coord_; size_t api_coord_n_; vector2ptr(coord, &api_coord_, &api_coord_n_);
        double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
        gmshModelMeshSetNode(nodeTag, api_coord_, api_coord_n_, api_parametricCoord_, api_parametricCoord_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_coord_);
        gmshFree(api_parametricCoord_);
      }

      // Rebuild the node cache.
      inline void rebuildNodeCache(const bool onlyIfNecessary = true)
      {
        int ierr = 0;
        gmshModelMeshRebuildNodeCache((int)onlyIfNecessary, &ierr);
        if(ierr) throwLastError();
      }

      // Rebuild the element cache.
      inline void rebuildElementCache(const bool onlyIfNecessary = true)
      {
        int ierr = 0;
        gmshModelMeshRebuildElementCache((int)onlyIfNecessary, &ierr);
        if(ierr) throwLastError();
      }

      // Get the nodes from all the elements belonging to the physical group of
      // dimension `dim' and tag `tag'. `nodeTags' contains the node tags; `coord'
      // is a vector of length 3 times the length of `nodeTags' that contains the
      // x, y, z coordinates of the nodes, concatenated: [n1x, n1y, n1z, n2x, ...].
      inline void getNodesForPhysicalGroup(const int dim,
                                           const int tag,
                                           std::vector<std::size_t> & nodeTags,
                                           std::vector<double> & coord)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        double *api_coord_; size_t api_coord_n_;
        gmshModelMeshGetNodesForPhysicalGroup(dim, tag, &api_nodeTags_, &api_nodeTags_n_, &api_coord_, &api_coord_n_, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
      }

      // Get the maximum tag `maxTag' of a node in the mesh.
      inline void getMaxNodeTag(std::size_t & maxTag)
      {
        int ierr = 0;
        gmshModelMeshGetMaxNodeTag(&maxTag, &ierr);
        if(ierr) throwLastError();
      }

      // Add nodes classified on the model entity of dimension `dim' and tag `tag'.
      // `nodeTags' contains the node tags (their unique, strictly positive
      // identification numbers). `coord' is a vector of length 3 times the length
      // of `nodeTags' that contains the x, y, z coordinates of the nodes,
      // concatenated: [n1x, n1y, n1z, n2x, ...]. The optional `parametricCoord'
      // vector contains the parametric coordinates of the nodes, if any. The
      // length of `parametricCoord' can be 0 or `dim' times the length of
      // `nodeTags'. If the `nodeTags' vector is empty, new tags are automatically
      // assigned to the nodes.
      inline void addNodes(const int dim,
                           const int tag,
                           const std::vector<std::size_t> & nodeTags,
                           const std::vector<double> & coord,
                           const std::vector<double> & parametricCoord = std::vector<double>())
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_; vector2ptr(nodeTags, &api_nodeTags_, &api_nodeTags_n_);
        double *api_coord_; size_t api_coord_n_; vector2ptr(coord, &api_coord_, &api_coord_n_);
        double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
        gmshModelMeshAddNodes(dim, tag, api_nodeTags_, api_nodeTags_n_, api_coord_, api_coord_n_, api_parametricCoord_, api_parametricCoord_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_nodeTags_);
        gmshFree(api_coord_);
        gmshFree(api_parametricCoord_);
      }

      // Reclassify all nodes on their associated model entity, based on the
      // elements. Can be used when importing nodes in bulk (e.g. by associating
      // them all to a single volume), to reclassify them correctly on model
      // surfaces, curves, etc. after the elements have been set.
      inline void reclassifyNodes()
      {
        int ierr = 0;
        gmshModelMeshReclassifyNodes(&ierr);
        if(ierr) throwLastError();
      }

      // Relocate the nodes classified on the entity of dimension `dim' and tag
      // `tag' using their parametric coordinates. If `tag' < 0, relocate the nodes
      // for all entities of dimension `dim'. If `dim' and `tag' are negative,
      // relocate all the nodes in the mesh.
      inline void relocateNodes(const int dim = -1,
                                const int tag = -1)
      {
        int ierr = 0;
        gmshModelMeshRelocateNodes(dim, tag, &ierr);
        if(ierr) throwLastError();
      }

      // Get the elements classified on the entity of dimension `dim' and tag
      // `tag'. If `tag' < 0, get the elements for all entities of dimension `dim'.
      // If `dim' and `tag' are negative, get all the elements in the mesh.
      // `elementTypes' contains the MSH types of the elements (e.g. `2' for 3-node
      // triangles: see `getElementProperties' to obtain the properties for a given
      // element type). `elementTags' is a vector of the same length as
      // `elementTypes'; each entry is a vector containing the tags (unique,
      // strictly positive identifiers) of the elements of the corresponding type.
      // `nodeTags' is also a vector of the same length as `elementTypes'; each
      // entry is a vector of length equal to the number of elements of the given
      // type times the number N of nodes for this type of element, that contains
      // the node tags of all the elements of the given type, concatenated: [e1n1,
      // e1n2, ..., e1nN, e2n1, ...].
      inline void getElements(std::vector<int> & elementTypes,
                              std::vector<std::vector<std::size_t> > & elementTags,
                              std::vector<std::vector<std::size_t> > & nodeTags,
                              const int dim = -1,
                              const int tag = -1)
      {
        int ierr = 0;
        int *api_elementTypes_; size_t api_elementTypes_n_;
        size_t **api_elementTags_; size_t *api_elementTags_n_, api_elementTags_nn_;
        size_t **api_nodeTags_; size_t *api_nodeTags_n_, api_nodeTags_nn_;
        gmshModelMeshGetElements(&api_elementTypes_, &api_elementTypes_n_, &api_elementTags_, &api_elementTags_n_, &api_elementTags_nn_, &api_nodeTags_, &api_nodeTags_n_, &api_nodeTags_nn_, dim, tag, &ierr);
        if(ierr) throwLastError();
        elementTypes.assign(api_elementTypes_, api_elementTypes_ + api_elementTypes_n_); gmshFree(api_elementTypes_);
        elementTags.resize(api_elementTags_nn_); for(size_t i = 0; i < api_elementTags_nn_; ++i){ elementTags[i].assign(api_elementTags_[i], api_elementTags_[i] + api_elementTags_n_[i]); gmshFree(api_elementTags_[i]); } gmshFree(api_elementTags_); gmshFree(api_elementTags_n_);
        nodeTags.resize(api_nodeTags_nn_); for(size_t i = 0; i < api_nodeTags_nn_; ++i){ nodeTags[i].assign(api_nodeTags_[i], api_nodeTags_[i] + api_nodeTags_n_[i]); gmshFree(api_nodeTags_[i]); } gmshFree(api_nodeTags_); gmshFree(api_nodeTags_n_);
      }

      // Get the type and node tags of the element with tag `tag', as well as the
      // dimension `dim' and tag `tag' of the entity on which the element is
      // classified. This function relies on an internal cache (a vector in case of
      // dense element numbering, a map otherwise); for large meshes accessing
      // elements in bulk is often preferable.
      inline void getElement(const std::size_t elementTag,
                             int & elementType,
                             std::vector<std::size_t> & nodeTags,
                             int & dim,
                             int & tag)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        gmshModelMeshGetElement(elementTag, &elementType, &api_nodeTags_, &api_nodeTags_n_, &dim, &tag, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
      }

      // Search the mesh for an element located at coordinates (`x', `y', `z').
      // This function performs a search in a spatial octree. If an element is
      // found, return its tag, type and node tags, as well as the local
      // coordinates (`u', `v', `w') within the reference element corresponding to
      // search location. If `dim' is >= 0, only search for elements of the given
      // dimension. If `strict' is not set, use a tolerance to find elements near
      // the search location.
      inline void getElementByCoordinates(const double x,
                                          const double y,
                                          const double z,
                                          std::size_t & elementTag,
                                          int & elementType,
                                          std::vector<std::size_t> & nodeTags,
                                          double & u,
                                          double & v,
                                          double & w,
                                          const int dim = -1,
                                          const bool strict = false)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        gmshModelMeshGetElementByCoordinates(x, y, z, &elementTag, &elementType, &api_nodeTags_, &api_nodeTags_n_, &u, &v, &w, dim, (int)strict, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
      }

      // Search the mesh for element(s) located at coordinates (`x', `y', `z').
      // This function performs a search in a spatial octree. Return the tags of
      // all found elements in `elementTags'. Additional information about the
      // elements can be accessed through `getElement' and
      // `getLocalCoordinatesInElement'. If `dim' is >= 0, only search for elements
      // of the given dimension. If `strict' is not set, use a tolerance to find
      // elements near the search location.
      inline void getElementsByCoordinates(const double x,
                                           const double y,
                                           const double z,
                                           std::vector<std::size_t> & elementTags,
                                           const int dim = -1,
                                           const bool strict = false)
      {
        int ierr = 0;
        size_t *api_elementTags_; size_t api_elementTags_n_;
        gmshModelMeshGetElementsByCoordinates(x, y, z, &api_elementTags_, &api_elementTags_n_, dim, (int)strict, &ierr);
        if(ierr) throwLastError();
        elementTags.assign(api_elementTags_, api_elementTags_ + api_elementTags_n_); gmshFree(api_elementTags_);
      }

      // Return the local coordinates (`u', `v', `w') within the element
      // `elementTag' corresponding to the model coordinates (`x', `y', `z'). This
      // function relies on an internal cache (a vector in case of dense element
      // numbering, a map otherwise); for large meshes accessing elements in bulk
      // is often preferable.
      inline void getLocalCoordinatesInElement(const std::size_t elementTag,
                                               const double x,
                                               const double y,
                                               const double z,
                                               double & u,
                                               double & v,
                                               double & w)
      {
        int ierr = 0;
        gmshModelMeshGetLocalCoordinatesInElement(elementTag, x, y, z, &u, &v, &w, &ierr);
        if(ierr) throwLastError();
      }

      // Get the types of elements in the entity of dimension `dim' and tag `tag'.
      // If `tag' < 0, get the types for all entities of dimension `dim'. If `dim'
      // and `tag' are negative, get all the types in the mesh.
      inline void getElementTypes(std::vector<int> & elementTypes,
                                  const int dim = -1,
                                  const int tag = -1)
      {
        int ierr = 0;
        int *api_elementTypes_; size_t api_elementTypes_n_;
        gmshModelMeshGetElementTypes(&api_elementTypes_, &api_elementTypes_n_, dim, tag, &ierr);
        if(ierr) throwLastError();
        elementTypes.assign(api_elementTypes_, api_elementTypes_ + api_elementTypes_n_); gmshFree(api_elementTypes_);
      }

      // Return an element type given its family name `familyName' ("Point",
      // "Line", "Triangle", "Quadrangle", "Tetrahedron", "Pyramid", "Prism",
      // "Hexahedron") and polynomial order `order'. If `serendip' is true, return
      // the corresponding serendip element type (element without interior nodes).
      inline int getElementType(const std::string & familyName,
                                const int order,
                                const bool serendip = false)
      {
        int ierr = 0;
        int result_api_ = gmshModelMeshGetElementType(familyName.c_str(), order, (int)serendip, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Get the properties of an element of type `elementType': its name
      // (`elementName'), dimension (`dim'), order (`order'), number of nodes
      // (`numNodes'), local coordinates of the nodes in the reference element
      // (`localNodeCoord' vector, of length `dim' times `numNodes') and number of
      // primary (first order) nodes (`numPrimaryNodes').
      inline void getElementProperties(const int elementType,
                                       std::string & elementName,
                                       int & dim,
                                       int & order,
                                       int & numNodes,
                                       std::vector<double> & localNodeCoord,
                                       int & numPrimaryNodes)
      {
        int ierr = 0;
        char *api_elementName_;
        double *api_localNodeCoord_; size_t api_localNodeCoord_n_;
        gmshModelMeshGetElementProperties(elementType, &api_elementName_, &dim, &order, &numNodes, &api_localNodeCoord_, &api_localNodeCoord_n_, &numPrimaryNodes, &ierr);
        if(ierr) throwLastError();
        elementName = std::string(api_elementName_); gmshFree(api_elementName_);
        localNodeCoord.assign(api_localNodeCoord_, api_localNodeCoord_ + api_localNodeCoord_n_); gmshFree(api_localNodeCoord_);
      }

      // Get the elements of type `elementType' classified on the entity of tag
      // `tag'. If `tag' < 0, get the elements for all entities. `elementTags' is a
      // vector containing the tags (unique, strictly positive identifiers) of the
      // elements of the corresponding type. `nodeTags' is a vector of length equal
      // to the number of elements of the given type times the number N of nodes
      // for this type of element, that contains the node tags of all the elements
      // of the given type, concatenated: [e1n1, e1n2, ..., e1nN, e2n1, ...]. If
      // `numTasks' > 1, only compute and return the part of the data indexed by
      // `task'.
      inline void getElementsByType(const int elementType,
                                    std::vector<std::size_t> & elementTags,
                                    std::vector<std::size_t> & nodeTags,
                                    const int tag = -1,
                                    const std::size_t task = 0,
                                    const std::size_t numTasks = 1)
      {
        int ierr = 0;
        size_t *api_elementTags_; size_t api_elementTags_n_;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        gmshModelMeshGetElementsByType(elementType, &api_elementTags_, &api_elementTags_n_, &api_nodeTags_, &api_nodeTags_n_, tag, task, numTasks, &ierr);
        if(ierr) throwLastError();
        elementTags.assign(api_elementTags_, api_elementTags_ + api_elementTags_n_); gmshFree(api_elementTags_);
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
      }

      // Get the maximum tag `maxTag' of an element in the mesh.
      inline void getMaxElementTag(std::size_t & maxTag)
      {
        int ierr = 0;
        gmshModelMeshGetMaxElementTag(&maxTag, &ierr);
        if(ierr) throwLastError();
      }

      // Preallocate data before calling `getElementsByType' with `numTasks' > 1.
      // For C and C++ only.
      inline void preallocateElementsByType(const int elementType,
                                            const bool elementTag,
                                            const bool nodeTag,
                                            std::vector<std::size_t> & elementTags,
                                            std::vector<std::size_t> & nodeTags,
                                            const int tag = -1)
      {
        int ierr = 0;
        size_t *api_elementTags_; size_t api_elementTags_n_;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        gmshModelMeshPreallocateElementsByType(elementType, (int)elementTag, (int)nodeTag, &api_elementTags_, &api_elementTags_n_, &api_nodeTags_, &api_nodeTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        elementTags.assign(api_elementTags_, api_elementTags_ + api_elementTags_n_); gmshFree(api_elementTags_);
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
      }

      // Add elements classified on the entity of dimension `dim' and tag `tag'.
      // `types' contains the MSH types of the elements (e.g. `2' for 3-node
      // triangles: see the Gmsh reference manual). `elementTags' is a vector of
      // the same length as `types'; each entry is a vector containing the tags
      // (unique, strictly positive identifiers) of the elements of the
      // corresponding type. `nodeTags' is also a vector of the same length as
      // `types'; each entry is a vector of length equal to the number of elements
      // of the given type times the number N of nodes per element, that contains
      // the node tags of all the elements of the given type, concatenated: [e1n1,
      // e1n2, ..., e1nN, e2n1, ...].
      inline void addElements(const int dim,
                              const int tag,
                              const std::vector<int> & elementTypes,
                              const std::vector<std::vector<std::size_t> > & elementTags,
                              const std::vector<std::vector<std::size_t> > & nodeTags)
      {
        int ierr = 0;
        int *api_elementTypes_; size_t api_elementTypes_n_; vector2ptr(elementTypes, &api_elementTypes_, &api_elementTypes_n_);
        size_t **api_elementTags_; size_t *api_elementTags_n_, api_elementTags_nn_; vectorvector2ptrptr(elementTags, &api_elementTags_, &api_elementTags_n_, &api_elementTags_nn_);
        size_t **api_nodeTags_; size_t *api_nodeTags_n_, api_nodeTags_nn_; vectorvector2ptrptr(nodeTags, &api_nodeTags_, &api_nodeTags_n_, &api_nodeTags_nn_);
        gmshModelMeshAddElements(dim, tag, api_elementTypes_, api_elementTypes_n_, (const size_t **)api_elementTags_, api_elementTags_n_, api_elementTags_nn_, (const size_t **)api_nodeTags_, api_nodeTags_n_, api_nodeTags_nn_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_elementTypes_);
        for(size_t i = 0; i < api_elementTags_nn_; ++i){ gmshFree(api_elementTags_[i]); } gmshFree(api_elementTags_); gmshFree(api_elementTags_n_);
        for(size_t i = 0; i < api_nodeTags_nn_; ++i){ gmshFree(api_nodeTags_[i]); } gmshFree(api_nodeTags_); gmshFree(api_nodeTags_n_);
      }

      // Add elements of type `elementType' classified on the entity of tag `tag'.
      // `elementTags' contains the tags (unique, strictly positive identifiers) of
      // the elements of the corresponding type. `nodeTags' is a vector of length
      // equal to the number of elements times the number N of nodes per element,
      // that contains the node tags of all the elements, concatenated: [e1n1,
      // e1n2, ..., e1nN, e2n1, ...]. If the `elementTag' vector is empty, new tags
      // are automatically assigned to the elements.
      inline void addElementsByType(const int tag,
                                    const int elementType,
                                    const std::vector<std::size_t> & elementTags,
                                    const std::vector<std::size_t> & nodeTags)
      {
        int ierr = 0;
        size_t *api_elementTags_; size_t api_elementTags_n_; vector2ptr(elementTags, &api_elementTags_, &api_elementTags_n_);
        size_t *api_nodeTags_; size_t api_nodeTags_n_; vector2ptr(nodeTags, &api_nodeTags_, &api_nodeTags_n_);
        gmshModelMeshAddElementsByType(tag, elementType, api_elementTags_, api_elementTags_n_, api_nodeTags_, api_nodeTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_elementTags_);
        gmshFree(api_nodeTags_);
      }

      // Get the numerical quadrature information for the given element type
      // `elementType' and integration rule `integrationType', where
      // `integrationType' concatenates the integration rule family name with the
      // desired order (e.g. "Gauss4" for a quadrature suited for integrating 4th
      // order polynomials). The "CompositeGauss" family uses tensor-product rules
      // based the 1D Gauss-Legendre rule; the "Gauss" family uses an economic
      // scheme when available (i.e. with a minimal number of points), and falls
      // back to "CompositeGauss" otherwise. Note that integration points for the
      // "Gauss" family can fall outside of the reference element for high-order
      // rules. `localCoord' contains the u, v, w coordinates of the G integration
      // points in the reference element: [g1u, g1v, g1w, ..., gGu, gGv, gGw].
      // `weights' contains the associated weights: [g1q, ..., gGq].
      inline void getIntegrationPoints(const int elementType,
                                       const std::string & integrationType,
                                       std::vector<double> & localCoord,
                                       std::vector<double> & weights)
      {
        int ierr = 0;
        double *api_localCoord_; size_t api_localCoord_n_;
        double *api_weights_; size_t api_weights_n_;
        gmshModelMeshGetIntegrationPoints(elementType, integrationType.c_str(), &api_localCoord_, &api_localCoord_n_, &api_weights_, &api_weights_n_, &ierr);
        if(ierr) throwLastError();
        localCoord.assign(api_localCoord_, api_localCoord_ + api_localCoord_n_); gmshFree(api_localCoord_);
        weights.assign(api_weights_, api_weights_ + api_weights_n_); gmshFree(api_weights_);
      }

      // Get the Jacobians of all the elements of type `elementType' classified on
      // the entity of tag `tag', at the G evaluation points `localCoord' given as
      // concatenated triplets of coordinates in the reference element [g1u, g1v,
      // g1w, ..., gGu, gGv, gGw]. Data is returned by element, with elements in
      // the same order as in `getElements' and `getElementsByType'. `jacobians'
      // contains for each element the 9 entries of the 3x3 Jacobian matrix at each
      // evaluation point. The matrix is returned by column: [e1g1Jxu, e1g1Jyu,
      // e1g1Jzu, e1g1Jxv, ..., e1g1Jzw, e1g2Jxu, ..., e1gGJzw, e2g1Jxu, ...], with
      // Jxu=dx/du, Jyu=dy/du, etc. `determinants' contains for each element the
      // determinant of the Jacobian matrix at each evaluation point: [e1g1, e1g2,
      // ... e1gG, e2g1, ...]. `coord' contains for each element the x, y, z
      // coordinates of the evaluation points. If `tag' < 0, get the Jacobian data
      // for all entities. If `numTasks' > 1, only compute and return the part of
      // the data indexed by `task'.
      inline void getJacobians(const int elementType,
                               const std::vector<double> & localCoord,
                               std::vector<double> & jacobians,
                               std::vector<double> & determinants,
                               std::vector<double> & coord,
                               const int tag = -1,
                               const std::size_t task = 0,
                               const std::size_t numTasks = 1)
      {
        int ierr = 0;
        double *api_localCoord_; size_t api_localCoord_n_; vector2ptr(localCoord, &api_localCoord_, &api_localCoord_n_);
        double *api_jacobians_; size_t api_jacobians_n_;
        double *api_determinants_; size_t api_determinants_n_;
        double *api_coord_; size_t api_coord_n_;
        gmshModelMeshGetJacobians(elementType, api_localCoord_, api_localCoord_n_, &api_jacobians_, &api_jacobians_n_, &api_determinants_, &api_determinants_n_, &api_coord_, &api_coord_n_, tag, task, numTasks, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_localCoord_);
        jacobians.assign(api_jacobians_, api_jacobians_ + api_jacobians_n_); gmshFree(api_jacobians_);
        determinants.assign(api_determinants_, api_determinants_ + api_determinants_n_); gmshFree(api_determinants_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
      }

      // Preallocate data before calling `getJacobians' with `numTasks' > 1. For C
      // and C++ only.
      inline void preallocateJacobians(const int elementType,
                                       const int numEvaluationPoints,
                                       const bool allocateJacobians,
                                       const bool allocateDeterminants,
                                       const bool allocateCoord,
                                       std::vector<double> & jacobians,
                                       std::vector<double> & determinants,
                                       std::vector<double> & coord,
                                       const int tag = -1)
      {
        int ierr = 0;
        double *api_jacobians_; size_t api_jacobians_n_;
        double *api_determinants_; size_t api_determinants_n_;
        double *api_coord_; size_t api_coord_n_;
        gmshModelMeshPreallocateJacobians(elementType, numEvaluationPoints, (int)allocateJacobians, (int)allocateDeterminants, (int)allocateCoord, &api_jacobians_, &api_jacobians_n_, &api_determinants_, &api_determinants_n_, &api_coord_, &api_coord_n_, tag, &ierr);
        if(ierr) throwLastError();
        jacobians.assign(api_jacobians_, api_jacobians_ + api_jacobians_n_); gmshFree(api_jacobians_);
        determinants.assign(api_determinants_, api_determinants_ + api_determinants_n_); gmshFree(api_determinants_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
      }

      // Get the Jacobian for a single element `elementTag', at the G evaluation
      // points `localCoord' given as concatenated triplets of coordinates in the
      // reference element [g1u, g1v, g1w, ..., gGu, gGv, gGw]. `jacobians'
      // contains the 9 entries of the 3x3 Jacobian matrix at each evaluation
      // point. The matrix is returned by column: [e1g1Jxu, e1g1Jyu, e1g1Jzu,
      // e1g1Jxv, ..., e1g1Jzw, e1g2Jxu, ..., e1gGJzw, e2g1Jxu, ...], with
      // Jxu=dx/du, Jyu=dy/du, etc. `determinants' contains the determinant of the
      // Jacobian matrix at each evaluation point. `coord' contains the x, y, z
      // coordinates of the evaluation points. This function relies on an internal
      // cache (a vector in case of dense element numbering, a map otherwise); for
      // large meshes accessing Jacobians in bulk is often preferable.
      inline void getJacobian(const std::size_t elementTag,
                              const std::vector<double> & localCoord,
                              std::vector<double> & jacobians,
                              std::vector<double> & determinants,
                              std::vector<double> & coord)
      {
        int ierr = 0;
        double *api_localCoord_; size_t api_localCoord_n_; vector2ptr(localCoord, &api_localCoord_, &api_localCoord_n_);
        double *api_jacobians_; size_t api_jacobians_n_;
        double *api_determinants_; size_t api_determinants_n_;
        double *api_coord_; size_t api_coord_n_;
        gmshModelMeshGetJacobian(elementTag, api_localCoord_, api_localCoord_n_, &api_jacobians_, &api_jacobians_n_, &api_determinants_, &api_determinants_n_, &api_coord_, &api_coord_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_localCoord_);
        jacobians.assign(api_jacobians_, api_jacobians_ + api_jacobians_n_); gmshFree(api_jacobians_);
        determinants.assign(api_determinants_, api_determinants_ + api_determinants_n_); gmshFree(api_determinants_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
      }

      // Get the basis functions of the element of type `elementType' at the
      // evaluation points `localCoord' (given as concatenated triplets of
      // coordinates in the reference element [g1u, g1v, g1w, ..., gGu, gGv, gGw]),
      // for the function space `functionSpaceType'. Currently supported function
      // spaces include "Lagrange" and "GradLagrange" for isoparametric Lagrange
      // basis functions and their gradient in the u, v, w coordinates of the
      // reference element; "LagrangeN" and "GradLagrangeN", with N = 1, 2, ...,
      // for N-th order Lagrange basis functions; "H1LegendreN" and
      // "GradH1LegendreN", with N = 1, 2, ..., for N-th order hierarchical H1
      // Legendre functions; "HcurlLegendreN" and "CurlHcurlLegendreN", with N = 1,
      // 2, ..., for N-th order curl-conforming basis functions. `numComponents'
      // returns the number C of components of a basis function (e.g. 1 for scalar
      // functions and 3 for vector functions). `basisFunctions' returns the value
      // of the N basis functions at the evaluation points, i.e. [g1f1, g1f2, ...,
      // g1fN, g2f1, ...] when C == 1 or [g1f1u, g1f1v, g1f1w, g1f2u, ..., g1fNw,
      // g2f1u, ...] when C == 3. For basis functions that depend on the
      // orientation of the elements, all values for the first orientation are
      // returned first, followed by values for the second, etc. `numOrientations'
      // returns the overall number of orientations. If `wantedOrientations' is not
      // empty, only return the values for the desired orientation indices.
      inline void getBasisFunctions(const int elementType,
                                    const std::vector<double> & localCoord,
                                    const std::string & functionSpaceType,
                                    int & numComponents,
                                    std::vector<double> & basisFunctions,
                                    int & numOrientations,
                                    const std::vector<int> & wantedOrientations = std::vector<int>())
      {
        int ierr = 0;
        double *api_localCoord_; size_t api_localCoord_n_; vector2ptr(localCoord, &api_localCoord_, &api_localCoord_n_);
        double *api_basisFunctions_; size_t api_basisFunctions_n_;
        int *api_wantedOrientations_; size_t api_wantedOrientations_n_; vector2ptr(wantedOrientations, &api_wantedOrientations_, &api_wantedOrientations_n_);
        gmshModelMeshGetBasisFunctions(elementType, api_localCoord_, api_localCoord_n_, functionSpaceType.c_str(), &numComponents, &api_basisFunctions_, &api_basisFunctions_n_, &numOrientations, api_wantedOrientations_, api_wantedOrientations_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_localCoord_);
        basisFunctions.assign(api_basisFunctions_, api_basisFunctions_ + api_basisFunctions_n_); gmshFree(api_basisFunctions_);
        gmshFree(api_wantedOrientations_);
      }

      // Get the orientation index of the elements of type `elementType' in the
      // entity of tag `tag'. The arguments have the same meaning as in
      // `getBasisFunctions'. `basisFunctionsOrientation' is a vector giving for
      // each element the orientation index in the values returned by
      // `getBasisFunctions'. For Lagrange basis functions the call is superfluous
      // as it will return a vector of zeros.
      inline void getBasisFunctionsOrientation(const int elementType,
                                               const std::string & functionSpaceType,
                                               std::vector<int> & basisFunctionsOrientation,
                                               const int tag = -1,
                                               const std::size_t task = 0,
                                               const std::size_t numTasks = 1)
      {
        int ierr = 0;
        int *api_basisFunctionsOrientation_; size_t api_basisFunctionsOrientation_n_;
        gmshModelMeshGetBasisFunctionsOrientation(elementType, functionSpaceType.c_str(), &api_basisFunctionsOrientation_, &api_basisFunctionsOrientation_n_, tag, task, numTasks, &ierr);
        if(ierr) throwLastError();
        basisFunctionsOrientation.assign(api_basisFunctionsOrientation_, api_basisFunctionsOrientation_ + api_basisFunctionsOrientation_n_); gmshFree(api_basisFunctionsOrientation_);
      }

      // Get the orientation of a single element `elementTag'.
      inline void getBasisFunctionsOrientationForElement(const std::size_t elementTag,
                                                         const std::string & functionSpaceType,
                                                         int & basisFunctionsOrientation)
      {
        int ierr = 0;
        gmshModelMeshGetBasisFunctionsOrientationForElement(elementTag, functionSpaceType.c_str(), &basisFunctionsOrientation, &ierr);
        if(ierr) throwLastError();
      }

      // Get the number of possible orientations for elements of type `elementType'
      // and function space named `functionSpaceType'.
      inline int getNumberOfOrientations(const int elementType,
                                         const std::string & functionSpaceType)
      {
        int ierr = 0;
        int result_api_ = gmshModelMeshGetNumberOfOrientations(elementType, functionSpaceType.c_str(), &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Preallocate data before calling `getBasisFunctionsOrientation' with
      // `numTasks' > 1. For C and C++ only.
      inline void preallocateBasisFunctionsOrientation(const int elementType,
                                                       std::vector<int> & basisFunctionsOrientation,
                                                       const int tag = -1)
      {
        int ierr = 0;
        int *api_basisFunctionsOrientation_; size_t api_basisFunctionsOrientation_n_;
        gmshModelMeshPreallocateBasisFunctionsOrientation(elementType, &api_basisFunctionsOrientation_, &api_basisFunctionsOrientation_n_, tag, &ierr);
        if(ierr) throwLastError();
        basisFunctionsOrientation.assign(api_basisFunctionsOrientation_, api_basisFunctionsOrientation_ + api_basisFunctionsOrientation_n_); gmshFree(api_basisFunctionsOrientation_);
      }

      // Get the global unique mesh edge identifiers `edgeTags' and orientations
      // `edgeOrientation' for an input list of node tag pairs defining these
      // edges, concatenated in the vector `nodeTags'. Mesh edges are created e.g.
      // by `createEdges()', `getKeys()' or `addEdges()'. The reference positive
      // orientation is n1 < n2, where n1 and n2 are the tags of the two edge
      // nodes, which corresponds to the local orientation of edge-based basis
      // functions as well.
      inline void getEdges(const std::vector<std::size_t> & nodeTags,
                           std::vector<std::size_t> & edgeTags,
                           std::vector<int> & edgeOrientations)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_; vector2ptr(nodeTags, &api_nodeTags_, &api_nodeTags_n_);
        size_t *api_edgeTags_; size_t api_edgeTags_n_;
        int *api_edgeOrientations_; size_t api_edgeOrientations_n_;
        gmshModelMeshGetEdges(api_nodeTags_, api_nodeTags_n_, &api_edgeTags_, &api_edgeTags_n_, &api_edgeOrientations_, &api_edgeOrientations_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_nodeTags_);
        edgeTags.assign(api_edgeTags_, api_edgeTags_ + api_edgeTags_n_); gmshFree(api_edgeTags_);
        edgeOrientations.assign(api_edgeOrientations_, api_edgeOrientations_ + api_edgeOrientations_n_); gmshFree(api_edgeOrientations_);
      }

      // Get the global unique mesh face identifiers `faceTags' and orientations
      // `faceOrientations' for an input list of node tag triplets (if `faceType'
      // == 3) or quadruplets (if `faceType' == 4) defining these faces,
      // concatenated in the vector `nodeTags'. Mesh faces are created e.g. by
      // `createFaces()', `getKeys()' or `addFaces()'.
      inline void getFaces(const int faceType,
                           const std::vector<std::size_t> & nodeTags,
                           std::vector<std::size_t> & faceTags,
                           std::vector<int> & faceOrientations)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_; vector2ptr(nodeTags, &api_nodeTags_, &api_nodeTags_n_);
        size_t *api_faceTags_; size_t api_faceTags_n_;
        int *api_faceOrientations_; size_t api_faceOrientations_n_;
        gmshModelMeshGetFaces(faceType, api_nodeTags_, api_nodeTags_n_, &api_faceTags_, &api_faceTags_n_, &api_faceOrientations_, &api_faceOrientations_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_nodeTags_);
        faceTags.assign(api_faceTags_, api_faceTags_ + api_faceTags_n_); gmshFree(api_faceTags_);
        faceOrientations.assign(api_faceOrientations_, api_faceOrientations_ + api_faceOrientations_n_); gmshFree(api_faceOrientations_);
      }

      // Create unique mesh edges for the entities `dimTags'.
      inline void createEdges(const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshCreateEdges(api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Create unique mesh faces for the entities `dimTags'.
      inline void createFaces(const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshCreateFaces(api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Get the global unique identifiers `edgeTags' and the nodes `edgeNodes' of
      // the edges in the mesh. Mesh edges are created e.g. by `createEdges()',
      // `getKeys()' or addEdges().
      inline void getAllEdges(std::vector<std::size_t> & edgeTags,
                              std::vector<std::size_t> & edgeNodes)
      {
        int ierr = 0;
        size_t *api_edgeTags_; size_t api_edgeTags_n_;
        size_t *api_edgeNodes_; size_t api_edgeNodes_n_;
        gmshModelMeshGetAllEdges(&api_edgeTags_, &api_edgeTags_n_, &api_edgeNodes_, &api_edgeNodes_n_, &ierr);
        if(ierr) throwLastError();
        edgeTags.assign(api_edgeTags_, api_edgeTags_ + api_edgeTags_n_); gmshFree(api_edgeTags_);
        edgeNodes.assign(api_edgeNodes_, api_edgeNodes_ + api_edgeNodes_n_); gmshFree(api_edgeNodes_);
      }

      // Get the global unique identifiers `faceTags' and the nodes `faceNodes' of
      // the faces of type `faceType' in the mesh. Mesh faces are created e.g. by
      // `createFaces()', `getKeys()' or addFaces().
      inline void getAllFaces(const int faceType,
                              std::vector<std::size_t> & faceTags,
                              std::vector<std::size_t> & faceNodes)
      {
        int ierr = 0;
        size_t *api_faceTags_; size_t api_faceTags_n_;
        size_t *api_faceNodes_; size_t api_faceNodes_n_;
        gmshModelMeshGetAllFaces(faceType, &api_faceTags_, &api_faceTags_n_, &api_faceNodes_, &api_faceNodes_n_, &ierr);
        if(ierr) throwLastError();
        faceTags.assign(api_faceTags_, api_faceTags_ + api_faceTags_n_); gmshFree(api_faceTags_);
        faceNodes.assign(api_faceNodes_, api_faceNodes_ + api_faceNodes_n_); gmshFree(api_faceNodes_);
      }

      // Add mesh edges defined by their global unique identifiers `edgeTags' and
      // their nodes `edgeNodes'.
      inline void addEdges(const std::vector<std::size_t> & edgeTags,
                           const std::vector<std::size_t> & edgeNodes)
      {
        int ierr = 0;
        size_t *api_edgeTags_; size_t api_edgeTags_n_; vector2ptr(edgeTags, &api_edgeTags_, &api_edgeTags_n_);
        size_t *api_edgeNodes_; size_t api_edgeNodes_n_; vector2ptr(edgeNodes, &api_edgeNodes_, &api_edgeNodes_n_);
        gmshModelMeshAddEdges(api_edgeTags_, api_edgeTags_n_, api_edgeNodes_, api_edgeNodes_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_edgeTags_);
        gmshFree(api_edgeNodes_);
      }

      // Add mesh faces of type `faceType' defined by their global unique
      // identifiers `faceTags' and their nodes `faceNodes'.
      inline void addFaces(const int faceType,
                           const std::vector<std::size_t> & faceTags,
                           const std::vector<std::size_t> & faceNodes)
      {
        int ierr = 0;
        size_t *api_faceTags_; size_t api_faceTags_n_; vector2ptr(faceTags, &api_faceTags_, &api_faceTags_n_);
        size_t *api_faceNodes_; size_t api_faceNodes_n_; vector2ptr(faceNodes, &api_faceNodes_, &api_faceNodes_n_);
        gmshModelMeshAddFaces(faceType, api_faceTags_, api_faceTags_n_, api_faceNodes_, api_faceNodes_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_faceTags_);
        gmshFree(api_faceNodes_);
      }

      // Generate the pair of keys for the elements of type `elementType' in the
      // entity of tag `tag', for the `functionSpaceType' function space. Each pair
      // (`typeKey', `entityKey') uniquely identifies a basis function in the
      // function space. If `returnCoord' is set, the `coord' vector contains the
      // x, y, z coordinates locating basis functions for sorting purposes.
      // Warning: this is an experimental feature and will probably change in a
      // future release.
      inline void getKeys(const int elementType,
                          const std::string & functionSpaceType,
                          std::vector<int> & typeKeys,
                          std::vector<std::size_t> & entityKeys,
                          std::vector<double> & coord,
                          const int tag = -1,
                          const bool returnCoord = true)
      {
        int ierr = 0;
        int *api_typeKeys_; size_t api_typeKeys_n_;
        size_t *api_entityKeys_; size_t api_entityKeys_n_;
        double *api_coord_; size_t api_coord_n_;
        gmshModelMeshGetKeys(elementType, functionSpaceType.c_str(), &api_typeKeys_, &api_typeKeys_n_, &api_entityKeys_, &api_entityKeys_n_, &api_coord_, &api_coord_n_, tag, (int)returnCoord, &ierr);
        if(ierr) throwLastError();
        typeKeys.assign(api_typeKeys_, api_typeKeys_ + api_typeKeys_n_); gmshFree(api_typeKeys_);
        entityKeys.assign(api_entityKeys_, api_entityKeys_ + api_entityKeys_n_); gmshFree(api_entityKeys_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
      }

      // Get the pair of keys for a single element `elementTag'.
      inline void getKeysForElement(const std::size_t elementTag,
                                    const std::string & functionSpaceType,
                                    std::vector<int> & typeKeys,
                                    std::vector<std::size_t> & entityKeys,
                                    std::vector<double> & coord,
                                    const bool returnCoord = true)
      {
        int ierr = 0;
        int *api_typeKeys_; size_t api_typeKeys_n_;
        size_t *api_entityKeys_; size_t api_entityKeys_n_;
        double *api_coord_; size_t api_coord_n_;
        gmshModelMeshGetKeysForElement(elementTag, functionSpaceType.c_str(), &api_typeKeys_, &api_typeKeys_n_, &api_entityKeys_, &api_entityKeys_n_, &api_coord_, &api_coord_n_, (int)returnCoord, &ierr);
        if(ierr) throwLastError();
        typeKeys.assign(api_typeKeys_, api_typeKeys_ + api_typeKeys_n_); gmshFree(api_typeKeys_);
        entityKeys.assign(api_entityKeys_, api_entityKeys_ + api_entityKeys_n_); gmshFree(api_entityKeys_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
      }

      // Get the number of keys by elements of type `elementType' for function
      // space named `functionSpaceType'.
      inline int getNumberOfKeys(const int elementType,
                                 const std::string & functionSpaceType)
      {
        int ierr = 0;
        int result_api_ = gmshModelMeshGetNumberOfKeys(elementType, functionSpaceType.c_str(), &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Get information about the pair of `keys'. `infoKeys' returns information
      // about the functions associated with the pairs (`typeKeys', `entityKey').
      // `infoKeys[0].first' describes the type of function (0 for  vertex
      // function, 1 for edge function, 2 for face function and 3 for bubble
      // function). `infoKeys[0].second' gives the order of the function associated
      // with the key. Warning: this is an experimental feature and will probably
      // change in a future release.
      inline void getKeysInformation(const std::vector<int> & typeKeys,
                                     const std::vector<std::size_t> & entityKeys,
                                     const int elementType,
                                     const std::string & functionSpaceType,
                                     gmsh::vectorpair & infoKeys)
      {
        int ierr = 0;
        int *api_typeKeys_; size_t api_typeKeys_n_; vector2ptr(typeKeys, &api_typeKeys_, &api_typeKeys_n_);
        size_t *api_entityKeys_; size_t api_entityKeys_n_; vector2ptr(entityKeys, &api_entityKeys_, &api_entityKeys_n_);
        int *api_infoKeys_; size_t api_infoKeys_n_;
        gmshModelMeshGetKeysInformation(api_typeKeys_, api_typeKeys_n_, api_entityKeys_, api_entityKeys_n_, elementType, functionSpaceType.c_str(), &api_infoKeys_, &api_infoKeys_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_typeKeys_);
        gmshFree(api_entityKeys_);
        infoKeys.resize(api_infoKeys_n_ / 2); for(size_t i = 0; i < api_infoKeys_n_ / 2; ++i){ infoKeys[i].first = api_infoKeys_[i * 2 + 0]; infoKeys[i].second = api_infoKeys_[i * 2 + 1]; } gmshFree(api_infoKeys_);
      }

      // Get the barycenters of all elements of type `elementType' classified on
      // the entity of tag `tag'. If `primary' is set, only the primary nodes of
      // the elements are taken into account for the barycenter calculation. If
      // `fast' is set, the function returns the sum of the primary node
      // coordinates (without normalizing by the number of nodes). If `tag' < 0,
      // get the barycenters for all entities. If `numTasks' > 1, only compute and
      // return the part of the data indexed by `task'.
      inline void getBarycenters(const int elementType,
                                 const int tag,
                                 const bool fast,
                                 const bool primary,
                                 std::vector<double> & barycenters,
                                 const std::size_t task = 0,
                                 const std::size_t numTasks = 1)
      {
        int ierr = 0;
        double *api_barycenters_; size_t api_barycenters_n_;
        gmshModelMeshGetBarycenters(elementType, tag, (int)fast, (int)primary, &api_barycenters_, &api_barycenters_n_, task, numTasks, &ierr);
        if(ierr) throwLastError();
        barycenters.assign(api_barycenters_, api_barycenters_ + api_barycenters_n_); gmshFree(api_barycenters_);
      }

      // Preallocate data before calling `getBarycenters' with `numTasks' > 1. For
      // C and C++ only.
      inline void preallocateBarycenters(const int elementType,
                                         std::vector<double> & barycenters,
                                         const int tag = -1)
      {
        int ierr = 0;
        double *api_barycenters_; size_t api_barycenters_n_;
        gmshModelMeshPreallocateBarycenters(elementType, &api_barycenters_, &api_barycenters_n_, tag, &ierr);
        if(ierr) throwLastError();
        barycenters.assign(api_barycenters_, api_barycenters_ + api_barycenters_n_); gmshFree(api_barycenters_);
      }

      // Get the nodes on the edges of all elements of type `elementType'
      // classified on the entity of tag `tag'. `nodeTags' contains the node tags
      // of the edges for all the elements: [e1a1n1, e1a1n2, e1a2n1, ...]. Data is
      // returned by element, with elements in the same order as in `getElements'
      // and `getElementsByType'. If `primary' is set, only the primary (begin/end)
      // nodes of the edges are returned. If `tag' < 0, get the edge nodes for all
      // entities. If `numTasks' > 1, only compute and return the part of the data
      // indexed by `task'.
      inline void getElementEdgeNodes(const int elementType,
                                      std::vector<std::size_t> & nodeTags,
                                      const int tag = -1,
                                      const bool primary = false,
                                      const std::size_t task = 0,
                                      const std::size_t numTasks = 1)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        gmshModelMeshGetElementEdgeNodes(elementType, &api_nodeTags_, &api_nodeTags_n_, tag, (int)primary, task, numTasks, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
      }

      // Get the nodes on the faces of type `faceType' (3 for triangular faces, 4
      // for quadrangular faces) of all elements of type `elementType' classified
      // on the entity of tag `tag'. `nodeTags' contains the node tags of the faces
      // for all elements: [e1f1n1, ..., e1f1nFaceType, e1f2n1, ...]. Data is
      // returned by element, with elements in the same order as in `getElements'
      // and `getElementsByType'. If `primary' is set, only the primary (corner)
      // nodes of the faces are returned. If `tag' < 0, get the face nodes for all
      // entities. If `numTasks' > 1, only compute and return the part of the data
      // indexed by `task'.
      inline void getElementFaceNodes(const int elementType,
                                      const int faceType,
                                      std::vector<std::size_t> & nodeTags,
                                      const int tag = -1,
                                      const bool primary = false,
                                      const std::size_t task = 0,
                                      const std::size_t numTasks = 1)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        gmshModelMeshGetElementFaceNodes(elementType, faceType, &api_nodeTags_, &api_nodeTags_n_, tag, (int)primary, task, numTasks, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
      }

      // Get the ghost elements `elementTags' and their associated `partitions'
      // stored in the ghost entity of dimension `dim' and tag `tag'.
      inline void getGhostElements(const int dim,
                                   const int tag,
                                   std::vector<std::size_t> & elementTags,
                                   std::vector<int> & partitions)
      {
        int ierr = 0;
        size_t *api_elementTags_; size_t api_elementTags_n_;
        int *api_partitions_; size_t api_partitions_n_;
        gmshModelMeshGetGhostElements(dim, tag, &api_elementTags_, &api_elementTags_n_, &api_partitions_, &api_partitions_n_, &ierr);
        if(ierr) throwLastError();
        elementTags.assign(api_elementTags_, api_elementTags_ + api_elementTags_n_); gmshFree(api_elementTags_);
        partitions.assign(api_partitions_, api_partitions_ + api_partitions_n_); gmshFree(api_partitions_);
      }

      // Set a mesh size constraint on the model entities `dimTags'. Currently only
      // entities of dimension 0 (points) are handled.
      inline void setSize(const gmsh::vectorpair & dimTags,
                          const double size)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshSetSize(api_dimTags_, api_dimTags_n_, size, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Get the mesh size constraints (if any) associated with the model entities
      // `dimTags'. A zero entry in the output `sizes' vector indicates that no
      // size constraint is specified on the corresponding entity.
      inline void getSizes(const gmsh::vectorpair & dimTags,
                           std::vector<double> & sizes)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        double *api_sizes_; size_t api_sizes_n_;
        gmshModelMeshGetSizes(api_dimTags_, api_dimTags_n_, &api_sizes_, &api_sizes_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        sizes.assign(api_sizes_, api_sizes_ + api_sizes_n_); gmshFree(api_sizes_);
      }

      // Set mesh size constraints at the given parametric points `parametricCoord'
      // on the model entity of dimension `dim' and tag `tag'. Currently only
      // entities of dimension 1 (lines) are handled.
      inline void setSizeAtParametricPoints(const int dim,
                                            const int tag,
                                            const std::vector<double> & parametricCoord,
                                            const std::vector<double> & sizes)
      {
        int ierr = 0;
        double *api_parametricCoord_; size_t api_parametricCoord_n_; vector2ptr(parametricCoord, &api_parametricCoord_, &api_parametricCoord_n_);
        double *api_sizes_; size_t api_sizes_n_; vector2ptr(sizes, &api_sizes_, &api_sizes_n_);
        gmshModelMeshSetSizeAtParametricPoints(dim, tag, api_parametricCoord_, api_parametricCoord_n_, api_sizes_, api_sizes_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_parametricCoord_);
        gmshFree(api_sizes_);
      }

      // Set a mesh size callback for the current model. The callback function
      // should take six arguments as input (`dim', `tag', `x', `y', `z' and `lc').
      // The first two integer arguments correspond to the dimension `dim' and tag
      // `tag' of the entity being meshed. The next four double precision arguments
      // correspond to the coordinates `x', `y' and `z' around which to prescribe
      // the mesh size and to the mesh size `lc' that would be prescribed if the
      // callback had not been called. The callback function should return a double
      // precision number specifying the desired mesh size; returning `lc' is
      // equivalent to a no-op.
      inline void setSizeCallback(std::function<double(int, int, double, double, double, double)> callback)
      {
        int ierr = 0;
        struct callback_caller_  {
          static double call(int dim, int tag, double x, double y, double z, double lc, void * callbackp_) {
            return (*static_cast<std::function<double(int, int, double, double, double, double)>*> (callbackp_))(dim, tag, x, y, z, lc);
          }
        };
        // FIXME memory leak
        auto *callback_ptr_ = new std::function<double(int, int, double, double, double, double)>(callback);
        gmshModelMeshSetSizeCallback(&callback_caller_::call, callback_ptr_, &ierr);
        if(ierr) throwLastError();
      }

      // Remove the mesh size callback from the current model.
      inline void removeSizeCallback()
      {
        int ierr = 0;
        gmshModelMeshRemoveSizeCallback(&ierr);
        if(ierr) throwLastError();
      }

      // Set a transfinite meshing constraint on the curve `tag', with `numNodes'
      // nodes distributed according to `meshType' and `coef'. Currently supported
      // types are "Progression" (geometrical progression with power `coef'),
      // "Bump" (refinement toward both extremities of the curve) and "Beta" (beta
      // law).
      inline void setTransfiniteCurve(const int tag,
                                      const int numNodes,
                                      const std::string & meshType = "Progression",
                                      const double coef = 1.)
      {
        int ierr = 0;
        gmshModelMeshSetTransfiniteCurve(tag, numNodes, meshType.c_str(), coef, &ierr);
        if(ierr) throwLastError();
      }

      // Set a transfinite meshing constraint on the surface `tag'. `arrangement'
      // describes the arrangement of the triangles when the surface is not flagged
      // as recombined: currently supported values are "Left", "Right",
      // "AlternateLeft" and "AlternateRight". `cornerTags' can be used to specify
      // the (3 or 4) corners of the transfinite interpolation explicitly;
      // specifying the corners explicitly is mandatory if the surface has more
      // that 3 or 4 points on its boundary.
      inline void setTransfiniteSurface(const int tag,
                                        const std::string & arrangement = "Left",
                                        const std::vector<int> & cornerTags = std::vector<int>())
      {
        int ierr = 0;
        int *api_cornerTags_; size_t api_cornerTags_n_; vector2ptr(cornerTags, &api_cornerTags_, &api_cornerTags_n_);
        gmshModelMeshSetTransfiniteSurface(tag, arrangement.c_str(), api_cornerTags_, api_cornerTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_cornerTags_);
      }

      // Set a transfinite meshing constraint on the surface `tag'. `cornerTags'
      // can be used to specify the (6 or 8) corners of the transfinite
      // interpolation explicitly.
      inline void setTransfiniteVolume(const int tag,
                                       const std::vector<int> & cornerTags = std::vector<int>())
      {
        int ierr = 0;
        int *api_cornerTags_; size_t api_cornerTags_n_; vector2ptr(cornerTags, &api_cornerTags_, &api_cornerTags_n_);
        gmshModelMeshSetTransfiniteVolume(tag, api_cornerTags_, api_cornerTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_cornerTags_);
      }

      // Set transfinite meshing constraints on the model entities in `dimTag'.
      // Transfinite meshing constraints are added to the curves of the
      // quadrangular surfaces and to the faces of 6-sided volumes. Quadragular
      // faces with a corner angle superior to `cornerAngle' (in radians) are
      // ignored. The number of points is automatically determined from the sizing
      // constraints. If `dimTag' is empty, the constraints are applied to all
      // entities in the model. If `recombine' is true, the recombine flag is
      // automatically set on the transfinite surfaces.
      inline void setTransfiniteAutomatic(const gmsh::vectorpair & dimTags = gmsh::vectorpair(),
                                          const double cornerAngle = 2.35,
                                          const bool recombine = true)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshSetTransfiniteAutomatic(api_dimTags_, api_dimTags_n_, cornerAngle, (int)recombine, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Set a recombination meshing constraint on the model entity of dimension
      // `dim' and tag `tag'. Currently only entities of dimension 2 (to recombine
      // triangles into quadrangles) are supported.
      inline void setRecombine(const int dim,
                               const int tag)
      {
        int ierr = 0;
        gmshModelMeshSetRecombine(dim, tag, &ierr);
        if(ierr) throwLastError();
      }

      // Set a smoothing meshing constraint on the model entity of dimension `dim'
      // and tag `tag'. `val' iterations of a Laplace smoother are applied.
      inline void setSmoothing(const int dim,
                               const int tag,
                               const int val)
      {
        int ierr = 0;
        gmshModelMeshSetSmoothing(dim, tag, val, &ierr);
        if(ierr) throwLastError();
      }

      // Set a reverse meshing constraint on the model entity of dimension `dim'
      // and tag `tag'. If `val' is true, the mesh orientation will be reversed
      // with respect to the natural mesh orientation (i.e. the orientation
      // consistent with the orientation of the geometry). If `val' is false, the
      // mesh is left as-is.
      inline void setReverse(const int dim,
                             const int tag,
                             const bool val = true)
      {
        int ierr = 0;
        gmshModelMeshSetReverse(dim, tag, (int)val, &ierr);
        if(ierr) throwLastError();
      }

      // Set the meshing algorithm on the model entity of dimension `dim' and tag
      // `tag'. Currently only supported for `dim' == 2.
      inline void setAlgorithm(const int dim,
                               const int tag,
                               const int val)
      {
        int ierr = 0;
        gmshModelMeshSetAlgorithm(dim, tag, val, &ierr);
        if(ierr) throwLastError();
      }

      // Force the mesh size to be extended from the boundary, or not, for the
      // model entity of dimension `dim' and tag `tag'. Currently only supported
      // for `dim' == 2.
      inline void setSizeFromBoundary(const int dim,
                                      const int tag,
                                      const int val)
      {
        int ierr = 0;
        gmshModelMeshSetSizeFromBoundary(dim, tag, val, &ierr);
        if(ierr) throwLastError();
      }

      // Set a compound meshing constraint on the model entities of dimension `dim'
      // and tags `tags'. During meshing, compound entities are treated as a single
      // discrete entity, which is automatically reparametrized.
      inline void setCompound(const int dim,
                              const std::vector<int> & tags)
      {
        int ierr = 0;
        int *api_tags_; size_t api_tags_n_; vector2ptr(tags, &api_tags_, &api_tags_n_);
        gmshModelMeshSetCompound(dim, api_tags_, api_tags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_tags_);
      }

      // Set meshing constraints on the bounding surfaces of the volume of tag
      // `tag' so that all surfaces are oriented with outward pointing normals; and
      // if a mesh already exists, reorient it. Currently only available with the
      // OpenCASCADE kernel, as it relies on the STL triangulation.
      inline void setOutwardOrientation(const int tag)
      {
        int ierr = 0;
        gmshModelMeshSetOutwardOrientation(tag, &ierr);
        if(ierr) throwLastError();
      }

      // Remove all meshing constraints from the model entities `dimTags'. If
      // `dimTags' is empty, remove all constraings.
      inline void removeConstraints(const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshRemoveConstraints(api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Embed the model entities of dimension `dim' and tags `tags' in the
      // (`inDim', `inTag') model entity. The dimension `dim' can 0, 1 or 2 and
      // must be strictly smaller than `inDim', which must be either 2 or 3. The
      // embedded entities should not intersect each other or be part of the
      // boundary of the entity `inTag', whose mesh will conform to the mesh of the
      // embedded entities. With the OpenCASCADE kernel, if the `fragment'
      // operation is applied to entities of different dimensions, the lower
      // dimensional entities will be automatically embedded in the higher
      // dimensional entities if they are not on their boundary.
      inline void embed(const int dim,
                        const std::vector<int> & tags,
                        const int inDim,
                        const int inTag)
      {
        int ierr = 0;
        int *api_tags_; size_t api_tags_n_; vector2ptr(tags, &api_tags_, &api_tags_n_);
        gmshModelMeshEmbed(dim, api_tags_, api_tags_n_, inDim, inTag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_tags_);
      }

      // Remove embedded entities from the model entities `dimTags'. if `dim' is >=
      // 0, only remove embedded entities of the given dimension (e.g. embedded
      // points if `dim' == 0).
      inline void removeEmbedded(const gmsh::vectorpair & dimTags,
                                 const int dim = -1)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshRemoveEmbedded(api_dimTags_, api_dimTags_n_, dim, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Get the entities (if any) embedded in the model entity of dimension `dim'
      // and tag `tag'.
      inline void getEmbedded(const int dim,
                              const int tag,
                              gmsh::vectorpair & dimTags)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_;
        gmshModelMeshGetEmbedded(dim, tag, &api_dimTags_, &api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        dimTags.resize(api_dimTags_n_ / 2); for(size_t i = 0; i < api_dimTags_n_ / 2; ++i){ dimTags[i].first = api_dimTags_[i * 2 + 0]; dimTags[i].second = api_dimTags_[i * 2 + 1]; } gmshFree(api_dimTags_);
      }

      // Reorder the elements of type `elementType' classified on the entity of tag
      // `tag' according to `ordering'.
      inline void reorderElements(const int elementType,
                                  const int tag,
                                  const std::vector<std::size_t> & ordering)
      {
        int ierr = 0;
        size_t *api_ordering_; size_t api_ordering_n_; vector2ptr(ordering, &api_ordering_, &api_ordering_n_);
        gmshModelMeshReorderElements(elementType, tag, api_ordering_, api_ordering_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_ordering_);
      }

      // Renumber the node tags in a continuous sequence.
      inline void renumberNodes()
      {
        int ierr = 0;
        gmshModelMeshRenumberNodes(&ierr);
        if(ierr) throwLastError();
      }

      // Renumber the element tags in a continuous sequence.
      inline void renumberElements()
      {
        int ierr = 0;
        gmshModelMeshRenumberElements(&ierr);
        if(ierr) throwLastError();
      }

      // Set the meshes of the entities of dimension `dim' and tag `tags' as
      // periodic copies of the meshes of entities `tagsMaster', using the affine
      // transformation specified in `affineTransformation' (16 entries of a 4x4
      // matrix, by row). If used after meshing, generate the periodic node
      // correspondence information assuming the meshes of entities `tags'
      // effectively match the meshes of entities `tagsMaster' (useful for
      // structured and extruded meshes). Currently only available for @code{dim}
      // == 1 and @code{dim} == 2.
      inline void setPeriodic(const int dim,
                              const std::vector<int> & tags,
                              const std::vector<int> & tagsMaster,
                              const std::vector<double> & affineTransform)
      {
        int ierr = 0;
        int *api_tags_; size_t api_tags_n_; vector2ptr(tags, &api_tags_, &api_tags_n_);
        int *api_tagsMaster_; size_t api_tagsMaster_n_; vector2ptr(tagsMaster, &api_tagsMaster_, &api_tagsMaster_n_);
        double *api_affineTransform_; size_t api_affineTransform_n_; vector2ptr(affineTransform, &api_affineTransform_, &api_affineTransform_n_);
        gmshModelMeshSetPeriodic(dim, api_tags_, api_tags_n_, api_tagsMaster_, api_tagsMaster_n_, api_affineTransform_, api_affineTransform_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_tags_);
        gmshFree(api_tagsMaster_);
        gmshFree(api_affineTransform_);
      }

      // Get master entities `tagsMaster' for the entities of dimension `dim' and
      // tags `tags'.
      inline void getPeriodic(const int dim,
                              const std::vector<int> & tags,
                              std::vector<int> & tagMaster)
      {
        int ierr = 0;
        int *api_tags_; size_t api_tags_n_; vector2ptr(tags, &api_tags_, &api_tags_n_);
        int *api_tagMaster_; size_t api_tagMaster_n_;
        gmshModelMeshGetPeriodic(dim, api_tags_, api_tags_n_, &api_tagMaster_, &api_tagMaster_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_tags_);
        tagMaster.assign(api_tagMaster_, api_tagMaster_ + api_tagMaster_n_); gmshFree(api_tagMaster_);
      }

      // Get the master entity `tagMaster', the node tags `nodeTags' and their
      // corresponding master node tags `nodeTagsMaster', and the affine transform
      // `affineTransform' for the entity of dimension `dim' and tag `tag'. If
      // `includeHighOrderNodes' is set, include high-order nodes in the returned
      // data.
      inline void getPeriodicNodes(const int dim,
                                   const int tag,
                                   int & tagMaster,
                                   std::vector<std::size_t> & nodeTags,
                                   std::vector<std::size_t> & nodeTagsMaster,
                                   std::vector<double> & affineTransform,
                                   const bool includeHighOrderNodes = false)
      {
        int ierr = 0;
        size_t *api_nodeTags_; size_t api_nodeTags_n_;
        size_t *api_nodeTagsMaster_; size_t api_nodeTagsMaster_n_;
        double *api_affineTransform_; size_t api_affineTransform_n_;
        gmshModelMeshGetPeriodicNodes(dim, tag, &tagMaster, &api_nodeTags_, &api_nodeTags_n_, &api_nodeTagsMaster_, &api_nodeTagsMaster_n_, &api_affineTransform_, &api_affineTransform_n_, (int)includeHighOrderNodes, &ierr);
        if(ierr) throwLastError();
        nodeTags.assign(api_nodeTags_, api_nodeTags_ + api_nodeTags_n_); gmshFree(api_nodeTags_);
        nodeTagsMaster.assign(api_nodeTagsMaster_, api_nodeTagsMaster_ + api_nodeTagsMaster_n_); gmshFree(api_nodeTagsMaster_);
        affineTransform.assign(api_affineTransform_, api_affineTransform_ + api_affineTransform_n_); gmshFree(api_affineTransform_);
      }

      // Get the master entity `tagMaster' and the key pairs (`typeKeyMaster',
      // `entityKeyMaster') corresponding to the entity `tag' and the key pairs
      // (`typeKey', `entityKey') for the elements of type `elementType' and
      // function space type `functionSpaceType'. If `returnCoord' is set, the
      // `coord' and `coordMaster' vectors contain the x, y, z coordinates locating
      // basis functions for sorting purposes.
      inline void getPeriodicKeys(const int elementType,
                                  const std::string & functionSpaceType,
                                  const int tag,
                                  int & tagMaster,
                                  std::vector<int> & typeKeys,
                                  std::vector<int> & typeKeysMaster,
                                  std::vector<std::size_t> & entityKeys,
                                  std::vector<std::size_t> & entityKeysMaster,
                                  std::vector<double> & coord,
                                  std::vector<double> & coordMaster,
                                  const bool returnCoord = true)
      {
        int ierr = 0;
        int *api_typeKeys_; size_t api_typeKeys_n_;
        int *api_typeKeysMaster_; size_t api_typeKeysMaster_n_;
        size_t *api_entityKeys_; size_t api_entityKeys_n_;
        size_t *api_entityKeysMaster_; size_t api_entityKeysMaster_n_;
        double *api_coord_; size_t api_coord_n_;
        double *api_coordMaster_; size_t api_coordMaster_n_;
        gmshModelMeshGetPeriodicKeys(elementType, functionSpaceType.c_str(), tag, &tagMaster, &api_typeKeys_, &api_typeKeys_n_, &api_typeKeysMaster_, &api_typeKeysMaster_n_, &api_entityKeys_, &api_entityKeys_n_, &api_entityKeysMaster_, &api_entityKeysMaster_n_, &api_coord_, &api_coord_n_, &api_coordMaster_, &api_coordMaster_n_, (int)returnCoord, &ierr);
        if(ierr) throwLastError();
        typeKeys.assign(api_typeKeys_, api_typeKeys_ + api_typeKeys_n_); gmshFree(api_typeKeys_);
        typeKeysMaster.assign(api_typeKeysMaster_, api_typeKeysMaster_ + api_typeKeysMaster_n_); gmshFree(api_typeKeysMaster_);
        entityKeys.assign(api_entityKeys_, api_entityKeys_ + api_entityKeys_n_); gmshFree(api_entityKeys_);
        entityKeysMaster.assign(api_entityKeysMaster_, api_entityKeysMaster_ + api_entityKeysMaster_n_); gmshFree(api_entityKeysMaster_);
        coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
        coordMaster.assign(api_coordMaster_, api_coordMaster_ + api_coordMaster_n_); gmshFree(api_coordMaster_);
      }

      // Remove duplicate nodes in the mesh of the current model.
      inline void removeDuplicateNodes()
      {
        int ierr = 0;
        gmshModelMeshRemoveDuplicateNodes(&ierr);
        if(ierr) throwLastError();
      }

      // Split (into two triangles) all quadrangles in surface `tag' whose quality
      // is lower than `quality'. If `tag' < 0, split quadrangles in all surfaces.
      inline void splitQuadrangles(const double quality = 1.,
                                   const int tag = -1)
      {
        int ierr = 0;
        gmshModelMeshSplitQuadrangles(quality, tag, &ierr);
        if(ierr) throwLastError();
      }

      // Classify ("color") the surface mesh based on the angle threshold `angle'
      // (in radians), and create new discrete surfaces, curves and points
      // accordingly. If `boundary' is set, also create discrete curves on the
      // boundary if the surface is open. If `forReparametrization' is set, create
      // curves and surfaces that can be reparametrized using a single map. If
      // `curveAngle' is less than Pi, also force curves to be split according to
      // `curveAngle'. If `exportDiscrete' is set, clear any built-in CAD kernel
      // entities and export the discrete entities in the built-in CAD kernel.
      inline void classifySurfaces(const double angle,
                                   const bool boundary = true,
                                   const bool forReparametrization = false,
                                   const double curveAngle = M_PI,
                                   const bool exportDiscrete = true)
      {
        int ierr = 0;
        gmshModelMeshClassifySurfaces(angle, (int)boundary, (int)forReparametrization, curveAngle, (int)exportDiscrete, &ierr);
        if(ierr) throwLastError();
      }

      // Create a geometry for the discrete entities `dimTags' (represented solely
      // by a mesh, without an underlying CAD description), i.e. create a
      // parametrization for discrete curves and surfaces, assuming that each can
      // be parametrized with a single map. If `dimTags' is empty, create a
      // geometry for all the discrete entities.
      inline void createGeometry(const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelMeshCreateGeometry(api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Create a boundary representation from the mesh if the model does not have
      // one (e.g. when imported from mesh file formats with no BRep representation
      // of the underlying model). If `makeSimplyConnected' is set, enforce simply
      // connected discrete surfaces and volumes. If `exportDiscrete' is set, clear
      // any built-in CAD kernel entities and export the discrete entities in the
      // built-in CAD kernel.
      inline void createTopology(const bool makeSimplyConnected = true,
                                 const bool exportDiscrete = true)
      {
        int ierr = 0;
        gmshModelMeshCreateTopology((int)makeSimplyConnected, (int)exportDiscrete, &ierr);
        if(ierr) throwLastError();
      }

      // Compute a basis representation for homology spaces after a mesh has been
      // generated. The computation domain is given in a list of physical group
      // tags `domainTags'; if empty, the whole mesh is the domain. The computation
      // subdomain for relative homology computation is given in a list of physical
      // group tags `subdomainTags'; if empty, absolute homology is computed. The
      // dimensions homology bases to be computed are given in the list `dim'; if
      // empty, all bases are computed. Resulting basis representation chains are
      // stored as physical groups in the mesh.
      inline void computeHomology(const std::vector<int> & domainTags = std::vector<int>(),
                                  const std::vector<int> & subdomainTags = std::vector<int>(),
                                  const std::vector<int> & dims = std::vector<int>())
      {
        int ierr = 0;
        int *api_domainTags_; size_t api_domainTags_n_; vector2ptr(domainTags, &api_domainTags_, &api_domainTags_n_);
        int *api_subdomainTags_; size_t api_subdomainTags_n_; vector2ptr(subdomainTags, &api_subdomainTags_, &api_subdomainTags_n_);
        int *api_dims_; size_t api_dims_n_; vector2ptr(dims, &api_dims_, &api_dims_n_);
        gmshModelMeshComputeHomology(api_domainTags_, api_domainTags_n_, api_subdomainTags_, api_subdomainTags_n_, api_dims_, api_dims_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_domainTags_);
        gmshFree(api_subdomainTags_);
        gmshFree(api_dims_);
      }

      // Compute a basis representation for cohomology spaces after a mesh has been
      // generated. The computation domain is given in a list of physical group
      // tags `domainTags'; if empty, the whole mesh is the domain. The computation
      // subdomain for relative cohomology computation is given in a list of
      // physical group tags `subdomainTags'; if empty, absolute cohomology is
      // computed. The dimensions homology bases to be computed are given in the
      // list `dim'; if empty, all bases are computed. Resulting basis
      // representation cochains are stored as physical groups in the mesh.
      inline void computeCohomology(const std::vector<int> & domainTags = std::vector<int>(),
                                    const std::vector<int> & subdomainTags = std::vector<int>(),
                                    const std::vector<int> & dims = std::vector<int>())
      {
        int ierr = 0;
        int *api_domainTags_; size_t api_domainTags_n_; vector2ptr(domainTags, &api_domainTags_, &api_domainTags_n_);
        int *api_subdomainTags_; size_t api_subdomainTags_n_; vector2ptr(subdomainTags, &api_subdomainTags_, &api_subdomainTags_n_);
        int *api_dims_; size_t api_dims_n_; vector2ptr(dims, &api_dims_, &api_dims_n_);
        gmshModelMeshComputeCohomology(api_domainTags_, api_domainTags_n_, api_subdomainTags_, api_subdomainTags_n_, api_dims_, api_dims_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_domainTags_);
        gmshFree(api_subdomainTags_);
        gmshFree(api_dims_);
      }

      // Compute a cross field for the current mesh. The function creates 3 views:
      // the H function, the Theta function and cross directions. Return the tags
      // of the views.
      inline void computeCrossField(std::vector<int> & viewTags)
      {
        int ierr = 0;
        int *api_viewTags_; size_t api_viewTags_n_;
        gmshModelMeshComputeCrossField(&api_viewTags_, &api_viewTags_n_, &ierr);
        if(ierr) throwLastError();
        viewTags.assign(api_viewTags_, api_viewTags_ + api_viewTags_n_); gmshFree(api_viewTags_);
      }

      // Triangulate the points given in the `coord' vector as pairs of u, v
      // coordinates, and return the node tags (with numbering starting at 1) of
      // the resulting triangles in `tri'.
      inline void triangulate(const std::vector<double> & coord,
                              std::vector<std::size_t> & tri)
      {
        int ierr = 0;
        double *api_coord_; size_t api_coord_n_; vector2ptr(coord, &api_coord_, &api_coord_n_);
        size_t *api_tri_; size_t api_tri_n_;
        gmshModelMeshTriangulate(api_coord_, api_coord_n_, &api_tri_, &api_tri_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_coord_);
        tri.assign(api_tri_, api_tri_ + api_tri_n_); gmshFree(api_tri_);
      }

      // Tetrahedralize the points given in the `coord' vector as triplets of x, y,
      // z coordinates, and return the node tags (with numbering starting at 1) of
      // the resulting tetrahedra in `tetra'.
      inline void tetrahedralize(const std::vector<double> & coord,
                                 std::vector<std::size_t> & tetra)
      {
        int ierr = 0;
        double *api_coord_; size_t api_coord_n_; vector2ptr(coord, &api_coord_, &api_coord_n_);
        size_t *api_tetra_; size_t api_tetra_n_;
        gmshModelMeshTetrahedralize(api_coord_, api_coord_n_, &api_tetra_, &api_tetra_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_coord_);
        tetra.assign(api_tetra_, api_tetra_ + api_tetra_n_); gmshFree(api_tetra_);
      }

      namespace field { // Mesh size field functions

        // Add a new mesh size field of type `fieldType'. If `tag' is positive,
        // assign the tag explicitly; otherwise a new tag is assigned
        // automatically. Return the field tag.
        inline int add(const std::string & fieldType,
                       const int tag = -1)
        {
          int ierr = 0;
          int result_api_ = gmshModelMeshFieldAdd(fieldType.c_str(), tag, &ierr);
          if(ierr) throwLastError();
          return result_api_;
        }

        // Remove the field with tag `tag'.
        inline void remove(const int tag)
        {
          int ierr = 0;
          gmshModelMeshFieldRemove(tag, &ierr);
          if(ierr) throwLastError();
        }

        // Get the list of all fields.
        inline void list(std::vector<int> & tags)
        {
          int ierr = 0;
          int *api_tags_; size_t api_tags_n_;
          gmshModelMeshFieldList(&api_tags_, &api_tags_n_, &ierr);
          if(ierr) throwLastError();
          tags.assign(api_tags_, api_tags_ + api_tags_n_); gmshFree(api_tags_);
        }

        // Get the type `fieldType' of the field with tag `tag'.
        inline void getType(const int tag,
                            std::string & fileType)
        {
          int ierr = 0;
          char *api_fileType_;
          gmshModelMeshFieldGetType(tag, &api_fileType_, &ierr);
          if(ierr) throwLastError();
          fileType = std::string(api_fileType_); gmshFree(api_fileType_);
        }

        // Set the numerical option `option' to value `value' for field `tag'.
        inline void setNumber(const int tag,
                              const std::string & option,
                              const double value)
        {
          int ierr = 0;
          gmshModelMeshFieldSetNumber(tag, option.c_str(), value, &ierr);
          if(ierr) throwLastError();
        }

        // Get the value of the numerical option `option' for field `tag'.
        inline void getNumber(const int tag,
                              const std::string & option,
                              double & value)
        {
          int ierr = 0;
          gmshModelMeshFieldGetNumber(tag, option.c_str(), &value, &ierr);
          if(ierr) throwLastError();
        }

        // Set the string option `option' to value `value' for field `tag'.
        inline void setString(const int tag,
                              const std::string & option,
                              const std::string & value)
        {
          int ierr = 0;
          gmshModelMeshFieldSetString(tag, option.c_str(), value.c_str(), &ierr);
          if(ierr) throwLastError();
        }

        // Get the value of the string option `option' for field `tag'.
        inline void getString(const int tag,
                              const std::string & option,
                              std::string & value)
        {
          int ierr = 0;
          char *api_value_;
          gmshModelMeshFieldGetString(tag, option.c_str(), &api_value_, &ierr);
          if(ierr) throwLastError();
          value = std::string(api_value_); gmshFree(api_value_);
        }

        // Set the numerical list option `option' to value `value' for field `tag'.
        inline void setNumbers(const int tag,
                               const std::string & option,
                               const std::vector<double> & value)
        {
          int ierr = 0;
          double *api_value_; size_t api_value_n_; vector2ptr(value, &api_value_, &api_value_n_);
          gmshModelMeshFieldSetNumbers(tag, option.c_str(), api_value_, api_value_n_, &ierr);
          if(ierr) throwLastError();
          gmshFree(api_value_);
        }

        // Get the value of the numerical list option `option' for field `tag'.
        inline void getNumbers(const int tag,
                               const std::string & option,
                               std::vector<double> & value)
        {
          int ierr = 0;
          double *api_value_; size_t api_value_n_;
          gmshModelMeshFieldGetNumbers(tag, option.c_str(), &api_value_, &api_value_n_, &ierr);
          if(ierr) throwLastError();
          value.assign(api_value_, api_value_ + api_value_n_); gmshFree(api_value_);
        }

        // Set the field `tag' as the background mesh size field.
        inline void setAsBackgroundMesh(const int tag)
        {
          int ierr = 0;
          gmshModelMeshFieldSetAsBackgroundMesh(tag, &ierr);
          if(ierr) throwLastError();
        }

        // Set the field `tag' as a boundary layer size field.
        inline void setAsBoundaryLayer(const int tag)
        {
          int ierr = 0;
          gmshModelMeshFieldSetAsBoundaryLayer(tag, &ierr);
          if(ierr) throwLastError();
        }

      } // namespace field

    } // namespace mesh

    namespace geo { // Built-in CAD kernel functions

      // Add a geometrical point in the built-in CAD representation, at coordinates
      // (`x', `y', `z'). If `meshSize' is > 0, add a meshing constraint at that
      // point. If `tag' is positive, set the tag explicitly; otherwise a new tag
      // is selected automatically. Return the tag of the point. (Note that the
      // point will be added in the current model only after `synchronize' is
      // called. This behavior holds for all the entities added in the geo module.)
      inline int addPoint(const double x,
                          const double y,
                          const double z,
                          const double meshSize = 0.,
                          const int tag = -1)
      {
        int ierr = 0;
        int result_api_ = gmshModelGeoAddPoint(x, y, z, meshSize, tag, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a straight line segment in the built-in CAD representation, between
      // the two points with tags `startTag' and `endTag'. If `tag' is positive,
      // set the tag explicitly; otherwise a new tag is selected automatically.
      // Return the tag of the line.
      inline int addLine(const int startTag,
                         const int endTag,
                         const int tag = -1)
      {
        int ierr = 0;
        int result_api_ = gmshModelGeoAddLine(startTag, endTag, tag, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a circle arc (strictly smaller than Pi) in the built-in CAD
      // representation, between the two points with tags `startTag' and `endTag',
      // and with center `centerTag'. If `tag' is positive, set the tag explicitly;
      // otherwise a new tag is selected automatically. If (`nx', `ny', `nz') !=
      // (0, 0, 0), explicitly set the plane of the circle arc. Return the tag of
      // the circle arc.
      inline int addCircleArc(const int startTag,
                              const int centerTag,
                              const int endTag,
                              const int tag = -1,
                              const double nx = 0.,
                              const double ny = 0.,
                              const double nz = 0.)
      {
        int ierr = 0;
        int result_api_ = gmshModelGeoAddCircleArc(startTag, centerTag, endTag, tag, nx, ny, nz, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add an ellipse arc (strictly smaller than Pi) in the built-in CAD
      // representation, between the two points `startTag' and `endTag', and with
      // center `centerTag' and major axis point `majorTag'. If `tag' is positive,
      // set the tag explicitly; otherwise a new tag is selected automatically. If
      // (`nx', `ny', `nz') != (0, 0, 0), explicitly set the plane of the circle
      // arc. Return the tag of the ellipse arc.
      inline int addEllipseArc(const int startTag,
                               const int centerTag,
                               const int majorTag,
                               const int endTag,
                               const int tag = -1,
                               const double nx = 0.,
                               const double ny = 0.,
                               const double nz = 0.)
      {
        int ierr = 0;
        int result_api_ = gmshModelGeoAddEllipseArc(startTag, centerTag, majorTag, endTag, tag, nx, ny, nz, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a spline (Catmull-Rom) curve in the built-in CAD representation, going
      // through the points `pointTags'. If `tag' is positive, set the tag
      // explicitly; otherwise a new tag is selected automatically. Create a
      // periodic curve if the first and last points are the same. Return the tag
      // of the spline curve.
      inline int addSpline(const std::vector<int> & pointTags,
                           const int tag = -1)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int result_api_ = gmshModelGeoAddSpline(api_pointTags_, api_pointTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        return result_api_;
      }

      // Add a cubic b-spline curve in the built-in CAD representation, with
      // `pointTags' control points. If `tag' is positive, set the tag explicitly;
      // otherwise a new tag is selected automatically. Creates a periodic curve if
      // the first and last points are the same. Return the tag of the b-spline
      // curve.
      inline int addBSpline(const std::vector<int> & pointTags,
                            const int tag = -1)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int result_api_ = gmshModelGeoAddBSpline(api_pointTags_, api_pointTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        return result_api_;
      }

      // Add a Bezier curve in the built-in CAD representation, with `pointTags'
      // control points. If `tag' is positive, set the tag explicitly; otherwise a
      // new tag is selected automatically.  Return the tag of the Bezier curve.
      inline int addBezier(const std::vector<int> & pointTags,
                           const int tag = -1)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int result_api_ = gmshModelGeoAddBezier(api_pointTags_, api_pointTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        return result_api_;
      }

      // Add a polyline curve in the built-in CAD representation, going through the
      // points `pointTags'. If `tag' is positive, set the tag explicitly;
      // otherwise a new tag is selected automatically. Create a periodic curve if
      // the first and last points are the same. Return the tag of the polyline
      // curve.
      inline int addPolyline(const std::vector<int> & pointTags,
                             const int tag = -1)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int result_api_ = gmshModelGeoAddPolyline(api_pointTags_, api_pointTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        return result_api_;
      }

      // Add a spline (Catmull-Rom) curve in the built-in CAD representation, going
      // through points sampling the curves in `curveTags'. The density of sampling
      // points on each curve is governed by `numIntervals'. If `tag' is positive,
      // set the tag explicitly; otherwise a new tag is selected automatically.
      // Return the tag of the spline.
      inline int addCompoundSpline(const std::vector<int> & curveTags,
                                   const int numIntervals = 5,
                                   const int tag = -1)
      {
        int ierr = 0;
        int *api_curveTags_; size_t api_curveTags_n_; vector2ptr(curveTags, &api_curveTags_, &api_curveTags_n_);
        int result_api_ = gmshModelGeoAddCompoundSpline(api_curveTags_, api_curveTags_n_, numIntervals, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_curveTags_);
        return result_api_;
      }

      // Add a b-spline curve in the built-in CAD representation, with control
      // points sampling the curves in `curveTags'. The density of sampling points
      // on each curve is governed by `numIntervals'. If `tag' is positive, set the
      // tag explicitly; otherwise a new tag is selected automatically. Return the
      // tag of the b-spline.
      inline int addCompoundBSpline(const std::vector<int> & curveTags,
                                    const int numIntervals = 20,
                                    const int tag = -1)
      {
        int ierr = 0;
        int *api_curveTags_; size_t api_curveTags_n_; vector2ptr(curveTags, &api_curveTags_, &api_curveTags_n_);
        int result_api_ = gmshModelGeoAddCompoundBSpline(api_curveTags_, api_curveTags_n_, numIntervals, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_curveTags_);
        return result_api_;
      }

      // Add a curve loop (a closed wire) in the built-in CAD representation,
      // formed by the curves `curveTags'. `curveTags' should contain (signed) tags
      // of model entities of dimension 1 forming a closed loop: a negative tag
      // signifies that the underlying curve is considered with reversed
      // orientation. If `tag' is positive, set the tag explicitly; otherwise a new
      // tag is selected automatically. If `reorient' is set, automatically
      // reorient the curves if necessary. Return the tag of the curve loop.
      inline int addCurveLoop(const std::vector<int> & curveTags,
                              const int tag = -1,
                              const bool reorient = false)
      {
        int ierr = 0;
        int *api_curveTags_; size_t api_curveTags_n_; vector2ptr(curveTags, &api_curveTags_, &api_curveTags_n_);
        int result_api_ = gmshModelGeoAddCurveLoop(api_curveTags_, api_curveTags_n_, tag, (int)reorient, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_curveTags_);
        return result_api_;
      }

      // Add curve loops in the built-in CAD representation based on the curves
      // `curveTags'. Return the `tags' of found curve loops, if any.
      inline void addCurveLoops(const std::vector<int> & curveTags,
                                std::vector<int> & tags)
      {
        int ierr = 0;
        int *api_curveTags_; size_t api_curveTags_n_; vector2ptr(curveTags, &api_curveTags_, &api_curveTags_n_);
        int *api_tags_; size_t api_tags_n_;
        gmshModelGeoAddCurveLoops(api_curveTags_, api_curveTags_n_, &api_tags_, &api_tags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_curveTags_);
        tags.assign(api_tags_, api_tags_ + api_tags_n_); gmshFree(api_tags_);
      }

      // Add a plane surface in the built-in CAD representation, defined by one or
      // more curve loops `wireTags'. The first curve loop defines the exterior
      // contour; additional curve loop define holes. If `tag' is positive, set the
      // tag explicitly; otherwise a new tag is selected automatically. Return the
      // tag of the surface.
      inline int addPlaneSurface(const std::vector<int> & wireTags,
                                 const int tag = -1)
      {
        int ierr = 0;
        int *api_wireTags_; size_t api_wireTags_n_; vector2ptr(wireTags, &api_wireTags_, &api_wireTags_n_);
        int result_api_ = gmshModelGeoAddPlaneSurface(api_wireTags_, api_wireTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_wireTags_);
        return result_api_;
      }

      // Add a surface in the built-in CAD representation, filling the curve loops
      // in `wireTags' using transfinite interpolation. Currently only a single
      // curve loop is supported; this curve loop should be composed by 3 or 4
      // curves only. If `tag' is positive, set the tag explicitly; otherwise a new
      // tag is selected automatically. Return the tag of the surface.
      inline int addSurfaceFilling(const std::vector<int> & wireTags,
                                   const int tag = -1,
                                   const int sphereCenterTag = -1)
      {
        int ierr = 0;
        int *api_wireTags_; size_t api_wireTags_n_; vector2ptr(wireTags, &api_wireTags_, &api_wireTags_n_);
        int result_api_ = gmshModelGeoAddSurfaceFilling(api_wireTags_, api_wireTags_n_, tag, sphereCenterTag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_wireTags_);
        return result_api_;
      }

      // Add a surface loop (a closed shell) formed by `surfaceTags' in the built-
      // in CAD representation.  If `tag' is positive, set the tag explicitly;
      // otherwise a new tag is selected automatically. Return the tag of the
      // shell.
      inline int addSurfaceLoop(const std::vector<int> & surfaceTags,
                                const int tag = -1)
      {
        int ierr = 0;
        int *api_surfaceTags_; size_t api_surfaceTags_n_; vector2ptr(surfaceTags, &api_surfaceTags_, &api_surfaceTags_n_);
        int result_api_ = gmshModelGeoAddSurfaceLoop(api_surfaceTags_, api_surfaceTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_surfaceTags_);
        return result_api_;
      }

      // Add a volume (a region) in the built-in CAD representation, defined by one
      // or more shells `shellTags'. The first surface loop defines the exterior
      // boundary; additional surface loop define holes. If `tag' is positive, set
      // the tag explicitly; otherwise a new tag is selected automatically. Return
      // the tag of the volume.
      inline int addVolume(const std::vector<int> & shellTags,
                           const int tag = -1)
      {
        int ierr = 0;
        int *api_shellTags_; size_t api_shellTags_n_; vector2ptr(shellTags, &api_shellTags_, &api_shellTags_n_);
        int result_api_ = gmshModelGeoAddVolume(api_shellTags_, api_shellTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_shellTags_);
        return result_api_;
      }

      // Extrude the entities `dimTags' in the built-in CAD representation, using a
      // translation along (`dx', `dy', `dz'). Return extruded entities in
      // `outDimTags'. If `numElements' is not empty, also extrude the mesh: the
      // entries in `numElements' give the number of elements in each layer. If
      // `height' is not empty, it provides the (cumulative) height of the
      // different layers, normalized to 1. If `recombine' is set, recombine the
      // mesh in the layers.
      inline void extrude(const gmsh::vectorpair & dimTags,
                          const double dx,
                          const double dy,
                          const double dz,
                          gmsh::vectorpair & outDimTags,
                          const std::vector<int> & numElements = std::vector<int>(),
                          const std::vector<double> & heights = std::vector<double>(),
                          const bool recombine = false)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int *api_numElements_; size_t api_numElements_n_; vector2ptr(numElements, &api_numElements_, &api_numElements_n_);
        double *api_heights_; size_t api_heights_n_; vector2ptr(heights, &api_heights_, &api_heights_n_);
        gmshModelGeoExtrude(api_dimTags_, api_dimTags_n_, dx, dy, dz, &api_outDimTags_, &api_outDimTags_n_, api_numElements_, api_numElements_n_, api_heights_, api_heights_n_, (int)recombine, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        gmshFree(api_numElements_);
        gmshFree(api_heights_);
      }

      // Extrude the entities `dimTags' in the built-in CAD representation, using a
      // rotation of `angle' radians around the axis of revolution defined by the
      // point (`x', `y', `z') and the direction (`ax', `ay', `az'). The angle
      // should be strictly smaller than Pi. Return extruded entities in
      // `outDimTags'. If `numElements' is not empty, also extrude the mesh: the
      // entries in `numElements' give the number of elements in each layer. If
      // `height' is not empty, it provides the (cumulative) height of the
      // different layers, normalized to 1. If `recombine' is set, recombine the
      // mesh in the layers.
      inline void revolve(const gmsh::vectorpair & dimTags,
                          const double x,
                          const double y,
                          const double z,
                          const double ax,
                          const double ay,
                          const double az,
                          const double angle,
                          gmsh::vectorpair & outDimTags,
                          const std::vector<int> & numElements = std::vector<int>(),
                          const std::vector<double> & heights = std::vector<double>(),
                          const bool recombine = false)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int *api_numElements_; size_t api_numElements_n_; vector2ptr(numElements, &api_numElements_, &api_numElements_n_);
        double *api_heights_; size_t api_heights_n_; vector2ptr(heights, &api_heights_, &api_heights_n_);
        gmshModelGeoRevolve(api_dimTags_, api_dimTags_n_, x, y, z, ax, ay, az, angle, &api_outDimTags_, &api_outDimTags_n_, api_numElements_, api_numElements_n_, api_heights_, api_heights_n_, (int)recombine, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        gmshFree(api_numElements_);
        gmshFree(api_heights_);
      }

      // Extrude the entities `dimTags' in the built-in CAD representation, using a
      // combined translation and rotation of `angle' radians, along (`dx', `dy',
      // `dz') and around the axis of revolution defined by the point (`x', `y',
      // `z') and the direction (`ax', `ay', `az'). The angle should be strictly
      // smaller than Pi. Return extruded entities in `outDimTags'. If
      // `numElements' is not empty, also extrude the mesh: the entries in
      // `numElements' give the number of elements in each layer. If `height' is
      // not empty, it provides the (cumulative) height of the different layers,
      // normalized to 1. If `recombine' is set, recombine the mesh in the layers.
      inline void twist(const gmsh::vectorpair & dimTags,
                        const double x,
                        const double y,
                        const double z,
                        const double dx,
                        const double dy,
                        const double dz,
                        const double ax,
                        const double ay,
                        const double az,
                        const double angle,
                        gmsh::vectorpair & outDimTags,
                        const std::vector<int> & numElements = std::vector<int>(),
                        const std::vector<double> & heights = std::vector<double>(),
                        const bool recombine = false)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int *api_numElements_; size_t api_numElements_n_; vector2ptr(numElements, &api_numElements_, &api_numElements_n_);
        double *api_heights_; size_t api_heights_n_; vector2ptr(heights, &api_heights_, &api_heights_n_);
        gmshModelGeoTwist(api_dimTags_, api_dimTags_n_, x, y, z, dx, dy, dz, ax, ay, az, angle, &api_outDimTags_, &api_outDimTags_n_, api_numElements_, api_numElements_n_, api_heights_, api_heights_n_, (int)recombine, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        gmshFree(api_numElements_);
        gmshFree(api_heights_);
      }

      // Extrude the entities `dimTags' in the built-in CAD representation along
      // the normals of the mesh, creating discrete boundary layer entities. Return
      // extruded entities in `outDimTags'. The entries in `numElements' give the
      // number of elements in each layer. If `height' is not empty, it provides
      // the height of the different layers. If `recombine' is set, recombine the
      // mesh in the layers. A second boundary layer can be created from the same
      // entities if `second' is set. If `viewIndex' is >= 0, use the corresponding
      // view to either specify the normals (if the view contains a vector field)
      // or scale the normals (if the view is scalar).
      inline void extrudeBoundaryLayer(const gmsh::vectorpair & dimTags,
                                       gmsh::vectorpair & outDimTags,
                                       const std::vector<int> & numElements = std::vector<int>(1, 1),
                                       const std::vector<double> & heights = std::vector<double>(),
                                       const bool recombine = false,
                                       const bool second = false,
                                       const int viewIndex = -1)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int *api_numElements_; size_t api_numElements_n_; vector2ptr(numElements, &api_numElements_, &api_numElements_n_);
        double *api_heights_; size_t api_heights_n_; vector2ptr(heights, &api_heights_, &api_heights_n_);
        gmshModelGeoExtrudeBoundaryLayer(api_dimTags_, api_dimTags_n_, &api_outDimTags_, &api_outDimTags_n_, api_numElements_, api_numElements_n_, api_heights_, api_heights_n_, (int)recombine, (int)second, viewIndex, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        gmshFree(api_numElements_);
        gmshFree(api_heights_);
      }

      // Translate the entities `dimTags' in the built-in CAD representation along
      // (`dx', `dy', `dz').
      inline void translate(const gmsh::vectorpair & dimTags,
                            const double dx,
                            const double dy,
                            const double dz)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelGeoTranslate(api_dimTags_, api_dimTags_n_, dx, dy, dz, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Rotate the entities `dimTags' in the built-in CAD representation by
      // `angle' radians around the axis of revolution defined by the point (`x',
      // `y', `z') and the direction (`ax', `ay', `az').
      inline void rotate(const gmsh::vectorpair & dimTags,
                         const double x,
                         const double y,
                         const double z,
                         const double ax,
                         const double ay,
                         const double az,
                         const double angle)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelGeoRotate(api_dimTags_, api_dimTags_n_, x, y, z, ax, ay, az, angle, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Scale the entities `dimTag' in the built-in CAD representation by factors
      // `a', `b' and `c' along the three coordinate axes; use (`x', `y', `z') as
      // the center of the homothetic transformation.
      inline void dilate(const gmsh::vectorpair & dimTags,
                         const double x,
                         const double y,
                         const double z,
                         const double a,
                         const double b,
                         const double c)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelGeoDilate(api_dimTags_, api_dimTags_n_, x, y, z, a, b, c, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Mirror the entities `dimTag' in the built-in CAD representation, with
      // respect to the plane of equation `a' * x + `b' * y + `c' * z + `d' = 0.
      inline void mirror(const gmsh::vectorpair & dimTags,
                         const double a,
                         const double b,
                         const double c,
                         const double d)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelGeoMirror(api_dimTags_, api_dimTags_n_, a, b, c, d, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Mirror the entities `dimTag' in the built-in CAD representation, with
      // respect to the plane of equation `a' * x + `b' * y + `c' * z + `d' = 0.
      // (This is a synonym for `mirror', which will be deprecated in a future
      // release.)
      inline void symmetrize(const gmsh::vectorpair & dimTags,
                             const double a,
                             const double b,
                             const double c,
                             const double d)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelGeoSymmetrize(api_dimTags_, api_dimTags_n_, a, b, c, d, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Copy the entities `dimTags' in the built-in CAD representation; the new
      // entities are returned in `outDimTags'.
      inline void copy(const gmsh::vectorpair & dimTags,
                       gmsh::vectorpair & outDimTags)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelGeoCopy(api_dimTags_, api_dimTags_n_, &api_outDimTags_, &api_outDimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Remove the entities `dimTags' in the built-in CAD representation, provided
      // that they are not on the boundary of higher-dimensional entities. If
      // `recursive' is true, remove all the entities on their boundaries, down to
      // dimension 0.
      inline void remove(const gmsh::vectorpair & dimTags,
                         const bool recursive = false)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelGeoRemove(api_dimTags_, api_dimTags_n_, (int)recursive, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Remove all duplicate entities in the built-in CAD representation
      // (different entities at the same geometrical location).
      inline void removeAllDuplicates()
      {
        int ierr = 0;
        gmshModelGeoRemoveAllDuplicates(&ierr);
        if(ierr) throwLastError();
      }

      // Split the curve of tag `tag' in the built-in CAD representation, on the
      // specified control points `pointTags'. This feature is only available for
      // lines, splines and b-splines. Return the tag(s) `curveTags' of the newly
      // created curve(s).
      inline void splitCurve(const int tag,
                             const std::vector<int> & pointTags,
                             std::vector<int> & curveTags)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int *api_curveTags_; size_t api_curveTags_n_;
        gmshModelGeoSplitCurve(tag, api_pointTags_, api_pointTags_n_, &api_curveTags_, &api_curveTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        curveTags.assign(api_curveTags_, api_curveTags_ + api_curveTags_n_); gmshFree(api_curveTags_);
      }

      // Get the maximum tag of entities of dimension `dim' in the built-in CAD
      // representation.
      inline int getMaxTag(const int dim)
      {
        int ierr = 0;
        int result_api_ = gmshModelGeoGetMaxTag(dim, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Set the maximum tag `maxTag' for entities of dimension `dim' in the built-
      // in CAD representation.
      inline void setMaxTag(const int dim,
                            const int maxTag)
      {
        int ierr = 0;
        gmshModelGeoSetMaxTag(dim, maxTag, &ierr);
        if(ierr) throwLastError();
      }

      // Add a physical group of dimension `dim', grouping the entities with tags
      // `tags' in the built-in CAD representation. Return the tag of the physical
      // group, equal to `tag' if `tag' is positive, or a new tag if `tag' < 0.
      inline int addPhysicalGroup(const int dim,
                                  const std::vector<int> & tags,
                                  const int tag = -1)
      {
        int ierr = 0;
        int *api_tags_; size_t api_tags_n_; vector2ptr(tags, &api_tags_, &api_tags_n_);
        int result_api_ = gmshModelGeoAddPhysicalGroup(dim, api_tags_, api_tags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_tags_);
        return result_api_;
      }

      // Remove the physical groups `dimTags' from the built-in CAD representation.
      // If `dimTags' is empty, remove all groups.
      inline void removePhysicalGroups(const gmsh::vectorpair & dimTags = gmsh::vectorpair())
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelGeoRemovePhysicalGroups(api_dimTags_, api_dimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Synchronize the built-in CAD representation with the current Gmsh model.
      // This can be called at any time, but since it involves a non trivial amount
      // of processing, the number of synchronization points should normally be
      // minimized. Without synchronization the entities in the built-in CAD
      // representation are not available to any function outside of the built-in
      // CAD kernel functions.
      inline void synchronize()
      {
        int ierr = 0;
        gmshModelGeoSynchronize(&ierr);
        if(ierr) throwLastError();
      }

      namespace mesh { // Built-in CAD kernel meshing constraints

        // Set a mesh size constraint on the entities `dimTags' in the built-in CAD
        // kernel representation. Currently only entities of dimension 0 (points)
        // are handled.
        inline void setSize(const gmsh::vectorpair & dimTags,
                            const double size)
        {
          int ierr = 0;
          int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
          gmshModelGeoMeshSetSize(api_dimTags_, api_dimTags_n_, size, &ierr);
          if(ierr) throwLastError();
          gmshFree(api_dimTags_);
        }

        // Set a transfinite meshing constraint on the curve `tag' in the built-in
        // CAD kernel representation, with `numNodes' nodes distributed according
        // to `meshType' and `coef'. Currently supported types are "Progression"
        // (geometrical progression with power `coef') and "Bump" (refinement
        // toward both extremities of the curve).
        inline void setTransfiniteCurve(const int tag,
                                        const int nPoints,
                                        const std::string & meshType = "Progression",
                                        const double coef = 1.)
        {
          int ierr = 0;
          gmshModelGeoMeshSetTransfiniteCurve(tag, nPoints, meshType.c_str(), coef, &ierr);
          if(ierr) throwLastError();
        }

        // Set a transfinite meshing constraint on the surface `tag' in the built-
        // in CAD kernel representation. `arrangement' describes the arrangement of
        // the triangles when the surface is not flagged as recombined: currently
        // supported values are "Left", "Right", "AlternateLeft" and
        // "AlternateRight". `cornerTags' can be used to specify the (3 or 4)
        // corners of the transfinite interpolation explicitly; specifying the
        // corners explicitly is mandatory if the surface has more that 3 or 4
        // points on its boundary.
        inline void setTransfiniteSurface(const int tag,
                                          const std::string & arrangement = "Left",
                                          const std::vector<int> & cornerTags = std::vector<int>())
        {
          int ierr = 0;
          int *api_cornerTags_; size_t api_cornerTags_n_; vector2ptr(cornerTags, &api_cornerTags_, &api_cornerTags_n_);
          gmshModelGeoMeshSetTransfiniteSurface(tag, arrangement.c_str(), api_cornerTags_, api_cornerTags_n_, &ierr);
          if(ierr) throwLastError();
          gmshFree(api_cornerTags_);
        }

        // Set a transfinite meshing constraint on the surface `tag' in the built-
        // in CAD kernel representation. `cornerTags' can be used to specify the (6
        // or 8) corners of the transfinite interpolation explicitly.
        inline void setTransfiniteVolume(const int tag,
                                         const std::vector<int> & cornerTags = std::vector<int>())
        {
          int ierr = 0;
          int *api_cornerTags_; size_t api_cornerTags_n_; vector2ptr(cornerTags, &api_cornerTags_, &api_cornerTags_n_);
          gmshModelGeoMeshSetTransfiniteVolume(tag, api_cornerTags_, api_cornerTags_n_, &ierr);
          if(ierr) throwLastError();
          gmshFree(api_cornerTags_);
        }

        // Set a recombination meshing constraint on the entity of dimension `dim'
        // and tag `tag' in the built-in CAD kernel representation. Currently only
        // entities of dimension 2 (to recombine triangles into quadrangles) are
        // supported.
        inline void setRecombine(const int dim,
                                 const int tag,
                                 const double angle = 45.)
        {
          int ierr = 0;
          gmshModelGeoMeshSetRecombine(dim, tag, angle, &ierr);
          if(ierr) throwLastError();
        }

        // Set a smoothing meshing constraint on the entity of dimension `dim' and
        // tag `tag' in the built-in CAD kernel representation. `val' iterations of
        // a Laplace smoother are applied.
        inline void setSmoothing(const int dim,
                                 const int tag,
                                 const int val)
        {
          int ierr = 0;
          gmshModelGeoMeshSetSmoothing(dim, tag, val, &ierr);
          if(ierr) throwLastError();
        }

        // Set a reverse meshing constraint on the entity of dimension `dim' and
        // tag `tag' in the built-in CAD kernel representation. If `val' is true,
        // the mesh orientation will be reversed with respect to the natural mesh
        // orientation (i.e. the orientation consistent with the orientation of the
        // geometry). If `val' is false, the mesh is left as-is.
        inline void setReverse(const int dim,
                               const int tag,
                               const bool val = true)
        {
          int ierr = 0;
          gmshModelGeoMeshSetReverse(dim, tag, (int)val, &ierr);
          if(ierr) throwLastError();
        }

        // Set the meshing algorithm on the entity of dimension `dim' and tag `tag'
        // in the built-in CAD kernel representation. Currently only supported for
        // `dim' == 2.
        inline void setAlgorithm(const int dim,
                                 const int tag,
                                 const int val)
        {
          int ierr = 0;
          gmshModelGeoMeshSetAlgorithm(dim, tag, val, &ierr);
          if(ierr) throwLastError();
        }

        // Force the mesh size to be extended from the boundary, or not, for the
        // entity of dimension `dim' and tag `tag' in the built-in CAD kernel
        // representation. Currently only supported for `dim' == 2.
        inline void setSizeFromBoundary(const int dim,
                                        const int tag,
                                        const int val)
        {
          int ierr = 0;
          gmshModelGeoMeshSetSizeFromBoundary(dim, tag, val, &ierr);
          if(ierr) throwLastError();
        }

      } // namespace mesh

    } // namespace geo

    namespace occ { // OpenCASCADE CAD kernel functions

      // Add a geometrical point in the OpenCASCADE CAD representation, at
      // coordinates (`x', `y', `z'). If `meshSize' is > 0, add a meshing
      // constraint at that point. If `tag' is positive, set the tag explicitly;
      // otherwise a new tag is selected automatically. Return the tag of the
      // point. (Note that the point will be added in the current model only after
      // `synchronize' is called. This behavior holds for all the entities added in
      // the occ module.)
      inline int addPoint(const double x,
                          const double y,
                          const double z,
                          const double meshSize = 0.,
                          const int tag = -1)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddPoint(x, y, z, meshSize, tag, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a straight line segment in the OpenCASCADE CAD representation, between
      // the two points with tags `startTag' and `endTag'. If `tag' is positive,
      // set the tag explicitly; otherwise a new tag is selected automatically.
      // Return the tag of the line.
      inline int addLine(const int startTag,
                         const int endTag,
                         const int tag = -1)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddLine(startTag, endTag, tag, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a circle arc in the OpenCASCADE CAD representation, between the two
      // points with tags `startTag' and `endTag', with center `centerTag'. If
      // `tag' is positive, set the tag explicitly; otherwise a new tag is selected
      // automatically. Return the tag of the circle arc.
      inline int addCircleArc(const int startTag,
                              const int centerTag,
                              const int endTag,
                              const int tag = -1)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddCircleArc(startTag, centerTag, endTag, tag, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a circle of center (`x', `y', `z') and radius `r' in the OpenCASCADE
      // CAD representation. If `tag' is positive, set the tag explicitly;
      // otherwise a new tag is selected automatically. If `angle1' and `angle2'
      // are specified, create a circle arc between the two angles. Return the tag
      // of the circle.
      inline int addCircle(const double x,
                           const double y,
                           const double z,
                           const double r,
                           const int tag = -1,
                           const double angle1 = 0.,
                           const double angle2 = 2*M_PI)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddCircle(x, y, z, r, tag, angle1, angle2, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add an ellipse arc in the OpenCASCADE CAD representation, between the two
      // points `startTag' and `endTag', and with center `centerTag' and major axis
      // point `majorTag'. If `tag' is positive, set the tag explicitly; otherwise
      // a new tag is selected automatically. Return the tag of the ellipse arc.
      // Note that OpenCASCADE does not allow creating ellipse arcs with the major
      // radius smaller than the minor radius.
      inline int addEllipseArc(const int startTag,
                               const int centerTag,
                               const int majorTag,
                               const int endTag,
                               const int tag = -1)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddEllipseArc(startTag, centerTag, majorTag, endTag, tag, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add an ellipse of center (`x', `y', `z') and radii `r1' and `r2' along the
      // x- and y-axes, respectively, in the OpenCASCADE CAD representation. If
      // `tag' is positive, set the tag explicitly; otherwise a new tag is selected
      // automatically. If `angle1' and `angle2' are specified, create an ellipse
      // arc between the two angles. Return the tag of the ellipse. Note that
      // OpenCASCADE does not allow creating ellipses with the major radius (along
      // the x-axis) smaller than or equal to the minor radius (along the y-axis):
      // rotate the shape or use `addCircle' in such cases.
      inline int addEllipse(const double x,
                            const double y,
                            const double z,
                            const double r1,
                            const double r2,
                            const int tag = -1,
                            const double angle1 = 0.,
                            const double angle2 = 2*M_PI)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddEllipse(x, y, z, r1, r2, tag, angle1, angle2, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a spline (C2 b-spline) curve in the OpenCASCADE CAD representation,
      // going through the points `pointTags'. If `tag' is positive, set the tag
      // explicitly; otherwise a new tag is selected automatically. Create a
      // periodic curve if the first and last points are the same. Return the tag
      // of the spline curve.
      inline int addSpline(const std::vector<int> & pointTags,
                           const int tag = -1)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int result_api_ = gmshModelOccAddSpline(api_pointTags_, api_pointTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        return result_api_;
      }

      // Add a b-spline curve of degree `degree' in the OpenCASCADE CAD
      // representation, with `pointTags' control points. If `weights', `knots' or
      // `multiplicities' are not provided, default parameters are computed
      // automatically. If `tag' is positive, set the tag explicitly; otherwise a
      // new tag is selected automatically. Create a periodic curve if the first
      // and last points are the same. Return the tag of the b-spline curve.
      inline int addBSpline(const std::vector<int> & pointTags,
                            const int tag = -1,
                            const int degree = 3,
                            const std::vector<double> & weights = std::vector<double>(),
                            const std::vector<double> & knots = std::vector<double>(),
                            const std::vector<int> & multiplicities = std::vector<int>())
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        double *api_weights_; size_t api_weights_n_; vector2ptr(weights, &api_weights_, &api_weights_n_);
        double *api_knots_; size_t api_knots_n_; vector2ptr(knots, &api_knots_, &api_knots_n_);
        int *api_multiplicities_; size_t api_multiplicities_n_; vector2ptr(multiplicities, &api_multiplicities_, &api_multiplicities_n_);
        int result_api_ = gmshModelOccAddBSpline(api_pointTags_, api_pointTags_n_, tag, degree, api_weights_, api_weights_n_, api_knots_, api_knots_n_, api_multiplicities_, api_multiplicities_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        gmshFree(api_weights_);
        gmshFree(api_knots_);
        gmshFree(api_multiplicities_);
        return result_api_;
      }

      // Add a Bezier curve in the OpenCASCADE CAD representation, with `pointTags'
      // control points. If `tag' is positive, set the tag explicitly; otherwise a
      // new tag is selected automatically. Return the tag of the Bezier curve.
      inline int addBezier(const std::vector<int> & pointTags,
                           const int tag = -1)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int result_api_ = gmshModelOccAddBezier(api_pointTags_, api_pointTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        return result_api_;
      }

      // Add a wire (open or closed) in the OpenCASCADE CAD representation, formed
      // by the curves `curveTags'. Note that an OpenCASCADE wire can be made of
      // curves that share geometrically identical (but topologically different)
      // points. If `tag' is positive, set the tag explicitly; otherwise a new tag
      // is selected automatically. Return the tag of the wire.
      inline int addWire(const std::vector<int> & curveTags,
                         const int tag = -1,
                         const bool checkClosed = false)
      {
        int ierr = 0;
        int *api_curveTags_; size_t api_curveTags_n_; vector2ptr(curveTags, &api_curveTags_, &api_curveTags_n_);
        int result_api_ = gmshModelOccAddWire(api_curveTags_, api_curveTags_n_, tag, (int)checkClosed, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_curveTags_);
        return result_api_;
      }

      // Add a curve loop (a closed wire) in the OpenCASCADE CAD representation,
      // formed by the curves `curveTags'. `curveTags' should contain tags of
      // curves forming a closed loop. Negative tags can be specified for
      // compatibility with the built-in kernel, but are simply ignored: the wire
      // is oriented according to the orientation of its first curve. Note that an
      // OpenCASCADE curve loop can be made of curves that share geometrically
      // identical (but topologically different) points. If `tag' is positive, set
      // the tag explicitly; otherwise a new tag is selected automatically. Return
      // the tag of the curve loop.
      inline int addCurveLoop(const std::vector<int> & curveTags,
                              const int tag = -1)
      {
        int ierr = 0;
        int *api_curveTags_; size_t api_curveTags_n_; vector2ptr(curveTags, &api_curveTags_, &api_curveTags_n_);
        int result_api_ = gmshModelOccAddCurveLoop(api_curveTags_, api_curveTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_curveTags_);
        return result_api_;
      }

      // Add a rectangle in the OpenCASCADE CAD representation, with lower left
      // corner at (`x', `y', `z') and upper right corner at (`x' + `dx', `y' +
      // `dy', `z'). If `tag' is positive, set the tag explicitly; otherwise a new
      // tag is selected automatically. Round the corners if `roundedRadius' is
      // nonzero. Return the tag of the rectangle.
      inline int addRectangle(const double x,
                              const double y,
                              const double z,
                              const double dx,
                              const double dy,
                              const int tag = -1,
                              const double roundedRadius = 0.)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddRectangle(x, y, z, dx, dy, tag, roundedRadius, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a disk in the OpenCASCADE CAD representation, with center (`xc', `yc',
      // `zc') and radius `rx' along the x-axis and `ry' along the y-axis. If `tag'
      // is positive, set the tag explicitly; otherwise a new tag is selected
      // automatically. Return the tag of the disk.
      inline int addDisk(const double xc,
                         const double yc,
                         const double zc,
                         const double rx,
                         const double ry,
                         const int tag = -1)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddDisk(xc, yc, zc, rx, ry, tag, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a plane surface in the OpenCASCADE CAD representation, defined by one
      // or more curve loops (or closed wires) `wireTags'. The first curve loop
      // defines the exterior contour; additional curve loop define holes. If `tag'
      // is positive, set the tag explicitly; otherwise a new tag is selected
      // automatically. Return the tag of the surface.
      inline int addPlaneSurface(const std::vector<int> & wireTags,
                                 const int tag = -1)
      {
        int ierr = 0;
        int *api_wireTags_; size_t api_wireTags_n_; vector2ptr(wireTags, &api_wireTags_, &api_wireTags_n_);
        int result_api_ = gmshModelOccAddPlaneSurface(api_wireTags_, api_wireTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_wireTags_);
        return result_api_;
      }

      // Add a surface in the OpenCASCADE CAD representation, filling the curve
      // loop `wireTag'. If `tag' is positive, set the tag explicitly; otherwise a
      // new tag is selected automatically. Return the tag of the surface. If
      // `pointTags' are provided, force the surface to pass through the given
      // points. The other optional arguments are `degree' (the degree of the
      // energy criterion to minimize for computing the deformation of the
      // surface), `numPointsOnCurves' (the average number of points for
      // discretisation of the bounding curves), `numIter' (the maximum number of
      // iterations of the optimization process), `anisotropic' (improve
      // performance when the ratio of the length along the two parametric
      // coordinates of the surface is high), `tol2d' (tolerance to the constraints
      // in the parametric plane of the surface), `tol3d' (the maximum distance
      // allowed between the support surface and the constraints), `tolAng' (the
      // maximum angle allowed between the normal of the surface and the
      // constraints), `tolCurv' (the maximum difference of curvature allowed
      // between the surface and the constraint), `maxDegree' (the highest degree
      // which the polynomial defining the filling surface can have) and,
      // `maxSegments' (the largest number of segments which the filling surface
      // can have).
      inline int addSurfaceFilling(const int wireTag,
                                   const int tag = -1,
                                   const std::vector<int> & pointTags = std::vector<int>(),
                                   const int degree = 3,
                                   const int numPointsOnCurves = 15,
                                   const int numIter = 2,
                                   const bool anisotropic = false,
                                   const double tol2d = 0.00001,
                                   const double tol3d = 0.0001,
                                   const double tolAng = 0.01,
                                   const double tolCurv = 0.1,
                                   const int maxDegree = 8,
                                   const int maxSegments = 9)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int result_api_ = gmshModelOccAddSurfaceFilling(wireTag, tag, api_pointTags_, api_pointTags_n_, degree, numPointsOnCurves, numIter, (int)anisotropic, tol2d, tol3d, tolAng, tolCurv, maxDegree, maxSegments, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        return result_api_;
      }

      // Add a BSpline surface in the OpenCASCADE CAD representation, filling the
      // curve loop `wireTag'. The curve loop should be made of 2, 3 or 4 BSpline
      // curves. The optional `type' argument specifies the type of filling:
      // "Stretch" creates the flattest patch, "Curved" (the default) creates the
      // most rounded patch, and "Coons" creates a rounded patch with less depth
      // than "Curved". If `tag' is positive, set the tag explicitly; otherwise a
      // new tag is selected automatically. Return the tag of the surface.
      inline int addBSplineFilling(const int wireTag,
                                   const int tag = -1,
                                   const std::string & type = "")
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddBSplineFilling(wireTag, tag, type.c_str(), &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a Bezier surface in the OpenCASCADE CAD representation, filling the
      // curve loop `wireTag'. The curve loop should be made of 2, 3 or 4 Bezier
      // curves. The optional `type' argument specifies the type of filling:
      // "Stretch" creates the flattest patch, "Curved" (the default) creates the
      // most rounded patch, and "Coons" creates a rounded patch with less depth
      // than "Curved". If `tag' is positive, set the tag explicitly; otherwise a
      // new tag is selected automatically. Return the tag of the surface.
      inline int addBezierFilling(const int wireTag,
                                  const int tag = -1,
                                  const std::string & type = "")
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddBezierFilling(wireTag, tag, type.c_str(), &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a b-spline surface of degree `degreeU' x `degreeV' in the OpenCASCADE
      // CAD representation, with `pointTags' control points given as a single
      // vector [Pu1v1, ... Pu`numPointsU'v1, Pu1v2, ...]. If `weights', `knotsU',
      // `knotsV', `multiplicitiesU' or `multiplicitiesV' are not provided, default
      // parameters are computed automatically. If `tag' is positive, set the tag
      // explicitly; otherwise a new tag is selected automatically. If `wireTags'
      // is provided, trim the b-spline patch using the provided wires: the first
      // wire defines the external contour, the others define holes. If `wire3D' is
      // set, consider wire curves as 3D curves and project them on the b-spline
      // surface; otherwise consider the wire curves as defined in the parametric
      // space of the surface. Return the tag of the b-spline surface.
      inline int addBSplineSurface(const std::vector<int> & pointTags,
                                   const int numPointsU,
                                   const int tag = -1,
                                   const int degreeU = 3,
                                   const int degreeV = 3,
                                   const std::vector<double> & weights = std::vector<double>(),
                                   const std::vector<double> & knotsU = std::vector<double>(),
                                   const std::vector<double> & knotsV = std::vector<double>(),
                                   const std::vector<int> & multiplicitiesU = std::vector<int>(),
                                   const std::vector<int> & multiplicitiesV = std::vector<int>(),
                                   const std::vector<int> & wireTags = std::vector<int>(),
                                   const bool wire3D = false)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        double *api_weights_; size_t api_weights_n_; vector2ptr(weights, &api_weights_, &api_weights_n_);
        double *api_knotsU_; size_t api_knotsU_n_; vector2ptr(knotsU, &api_knotsU_, &api_knotsU_n_);
        double *api_knotsV_; size_t api_knotsV_n_; vector2ptr(knotsV, &api_knotsV_, &api_knotsV_n_);
        int *api_multiplicitiesU_; size_t api_multiplicitiesU_n_; vector2ptr(multiplicitiesU, &api_multiplicitiesU_, &api_multiplicitiesU_n_);
        int *api_multiplicitiesV_; size_t api_multiplicitiesV_n_; vector2ptr(multiplicitiesV, &api_multiplicitiesV_, &api_multiplicitiesV_n_);
        int *api_wireTags_; size_t api_wireTags_n_; vector2ptr(wireTags, &api_wireTags_, &api_wireTags_n_);
        int result_api_ = gmshModelOccAddBSplineSurface(api_pointTags_, api_pointTags_n_, numPointsU, tag, degreeU, degreeV, api_weights_, api_weights_n_, api_knotsU_, api_knotsU_n_, api_knotsV_, api_knotsV_n_, api_multiplicitiesU_, api_multiplicitiesU_n_, api_multiplicitiesV_, api_multiplicitiesV_n_, api_wireTags_, api_wireTags_n_, (int)wire3D, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        gmshFree(api_weights_);
        gmshFree(api_knotsU_);
        gmshFree(api_knotsV_);
        gmshFree(api_multiplicitiesU_);
        gmshFree(api_multiplicitiesV_);
        gmshFree(api_wireTags_);
        return result_api_;
      }

      // Add a Bezier surface in the OpenCASCADE CAD representation, with
      // `pointTags' control points given as a single vector [Pu1v1, ...
      // Pu`numPointsU'v1, Pu1v2, ...]. If `tag' is positive, set the tag
      // explicitly; otherwise a new tag is selected automatically. If `wireTags'
      // is provided, trim the Bezier patch using the provided wires: the first
      // wire defines the external contour, the others define holes. If `wire3D' is
      // set, consider wire curves as 3D curves and project them on the Bezier
      // surface; otherwise consider the wire curves as defined in the parametric
      // space of the surface. Return the tag of the Bezier surface.
      inline int addBezierSurface(const std::vector<int> & pointTags,
                                  const int numPointsU,
                                  const int tag = -1,
                                  const std::vector<int> & wireTags = std::vector<int>(),
                                  const bool wire3D = false)
      {
        int ierr = 0;
        int *api_pointTags_; size_t api_pointTags_n_; vector2ptr(pointTags, &api_pointTags_, &api_pointTags_n_);
        int *api_wireTags_; size_t api_wireTags_n_; vector2ptr(wireTags, &api_wireTags_, &api_wireTags_n_);
        int result_api_ = gmshModelOccAddBezierSurface(api_pointTags_, api_pointTags_n_, numPointsU, tag, api_wireTags_, api_wireTags_n_, (int)wire3D, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_pointTags_);
        gmshFree(api_wireTags_);
        return result_api_;
      }

      // Trim the surface `surfaceTag' with the wires `wireTags', replacing any
      // existing trimming curves. The first wire defines the external contour, the
      // others define holes. If `wire3D' is set, consider wire curves as 3D curves
      // and project them on the surface; otherwise consider the wire curves as
      // defined in the parametric space of the surface. If `tag' is positive, set
      // the tag explicitly; otherwise a new tag is selected automatically. Return
      // the tag of the trimmed surface.
      inline int addTrimmedSurface(const int surfaceTag,
                                   const std::vector<int> & wireTags = std::vector<int>(),
                                   const bool wire3D = false,
                                   const int tag = -1)
      {
        int ierr = 0;
        int *api_wireTags_; size_t api_wireTags_n_; vector2ptr(wireTags, &api_wireTags_, &api_wireTags_n_);
        int result_api_ = gmshModelOccAddTrimmedSurface(surfaceTag, api_wireTags_, api_wireTags_n_, (int)wire3D, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_wireTags_);
        return result_api_;
      }

      // Add a surface loop (a closed shell) in the OpenCASCADE CAD representation,
      // formed by `surfaceTags'.  If `tag' is positive, set the tag explicitly;
      // otherwise a new tag is selected automatically. Return the tag of the
      // surface loop. Setting `sewing' allows one to build a shell made of
      // surfaces that share geometrically identical (but topologically different)
      // curves.
      inline int addSurfaceLoop(const std::vector<int> & surfaceTags,
                                const int tag = -1,
                                const bool sewing = false)
      {
        int ierr = 0;
        int *api_surfaceTags_; size_t api_surfaceTags_n_; vector2ptr(surfaceTags, &api_surfaceTags_, &api_surfaceTags_n_);
        int result_api_ = gmshModelOccAddSurfaceLoop(api_surfaceTags_, api_surfaceTags_n_, tag, (int)sewing, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_surfaceTags_);
        return result_api_;
      }

      // Add a volume (a region) in the OpenCASCADE CAD representation, defined by
      // one or more surface loops `shellTags'. The first surface loop defines the
      // exterior boundary; additional surface loop define holes. If `tag' is
      // positive, set the tag explicitly; otherwise a new tag is selected
      // automatically. Return the tag of the volume.
      inline int addVolume(const std::vector<int> & shellTags,
                           const int tag = -1)
      {
        int ierr = 0;
        int *api_shellTags_; size_t api_shellTags_n_; vector2ptr(shellTags, &api_shellTags_, &api_shellTags_n_);
        int result_api_ = gmshModelOccAddVolume(api_shellTags_, api_shellTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_shellTags_);
        return result_api_;
      }

      // Add a sphere of center (`xc', `yc', `zc') and radius `r' in the
      // OpenCASCADE CAD representation. The optional `angle1' and `angle2'
      // arguments define the polar angle opening (from -Pi/2 to Pi/2). The
      // optional `angle3' argument defines the azimuthal opening (from 0 to 2*Pi).
      // If `tag' is positive, set the tag explicitly; otherwise a new tag is
      // selected automatically. Return the tag of the sphere.
      inline int addSphere(const double xc,
                           const double yc,
                           const double zc,
                           const double radius,
                           const int tag = -1,
                           const double angle1 = -M_PI/2,
                           const double angle2 = M_PI/2,
                           const double angle3 = 2*M_PI)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddSphere(xc, yc, zc, radius, tag, angle1, angle2, angle3, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a parallelepipedic box in the OpenCASCADE CAD representation, defined
      // by a point (`x', `y', `z') and the extents along the x-, y- and z-axes. If
      // `tag' is positive, set the tag explicitly; otherwise a new tag is selected
      // automatically. Return the tag of the box.
      inline int addBox(const double x,
                        const double y,
                        const double z,
                        const double dx,
                        const double dy,
                        const double dz,
                        const int tag = -1)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddBox(x, y, z, dx, dy, dz, tag, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a cylinder in the OpenCASCADE CAD representation, defined by the
      // center (`x', `y', `z') of its first circular face, the 3 components (`dx',
      // `dy', `dz') of the vector defining its axis and its radius `r'. The
      // optional `angle' argument defines the angular opening (from 0 to 2*Pi). If
      // `tag' is positive, set the tag explicitly; otherwise a new tag is selected
      // automatically. Return the tag of the cylinder.
      inline int addCylinder(const double x,
                             const double y,
                             const double z,
                             const double dx,
                             const double dy,
                             const double dz,
                             const double r,
                             const int tag = -1,
                             const double angle = 2*M_PI)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddCylinder(x, y, z, dx, dy, dz, r, tag, angle, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a cone in the OpenCASCADE CAD representation, defined by the center
      // (`x', `y', `z') of its first circular face, the 3 components of the vector
      // (`dx', `dy', `dz') defining its axis and the two radii `r1' and `r2' of
      // the faces (these radii can be zero). If `tag' is positive, set the tag
      // explicitly; otherwise a new tag is selected automatically. `angle' defines
      // the optional angular opening (from 0 to 2*Pi). Return the tag of the cone.
      inline int addCone(const double x,
                         const double y,
                         const double z,
                         const double dx,
                         const double dy,
                         const double dz,
                         const double r1,
                         const double r2,
                         const int tag = -1,
                         const double angle = 2*M_PI)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddCone(x, y, z, dx, dy, dz, r1, r2, tag, angle, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a right angular wedge in the OpenCASCADE CAD representation, defined
      // by the right-angle point (`x', `y', `z') and the 3 extends along the x-,
      // y- and z-axes (`dx', `dy', `dz'). If `tag' is positive, set the tag
      // explicitly; otherwise a new tag is selected automatically. The optional
      // argument `ltx' defines the top extent along the x-axis. Return the tag of
      // the wedge.
      inline int addWedge(const double x,
                          const double y,
                          const double z,
                          const double dx,
                          const double dy,
                          const double dz,
                          const int tag = -1,
                          const double ltx = 0.)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddWedge(x, y, z, dx, dy, dz, tag, ltx, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a torus in the OpenCASCADE CAD representation, defined by its center
      // (`x', `y', `z') and its 2 radii `r' and `r2'. If `tag' is positive, set
      // the tag explicitly; otherwise a new tag is selected automatically. The
      // optional argument `angle' defines the angular opening (from 0 to 2*Pi).
      // Return the tag of the wedge.
      inline int addTorus(const double x,
                          const double y,
                          const double z,
                          const double r1,
                          const double r2,
                          const int tag = -1,
                          const double angle = 2*M_PI)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccAddTorus(x, y, z, r1, r2, tag, angle, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Add a volume (if the optional argument `makeSolid' is set) or surfaces in
      // the OpenCASCADE CAD representation, defined through the open or closed
      // wires `wireTags'. If `tag' is positive, set the tag explicitly; otherwise
      // a new tag is selected automatically. The new entities are returned in
      // `outDimTags'. If the optional argument `makeRuled' is set, the surfaces
      // created on the boundary are forced to be ruled surfaces. If `maxDegree' is
      // positive, set the maximal degree of resulting surface.
      inline void addThruSections(const std::vector<int> & wireTags,
                                  gmsh::vectorpair & outDimTags,
                                  const int tag = -1,
                                  const bool makeSolid = true,
                                  const bool makeRuled = false,
                                  const int maxDegree = -1)
      {
        int ierr = 0;
        int *api_wireTags_; size_t api_wireTags_n_; vector2ptr(wireTags, &api_wireTags_, &api_wireTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelOccAddThruSections(api_wireTags_, api_wireTags_n_, &api_outDimTags_, &api_outDimTags_n_, tag, (int)makeSolid, (int)makeRuled, maxDegree, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_wireTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Add a hollowed volume in the OpenCASCADE CAD representation, built from an
      // initial volume `volumeTag' and a set of faces from this volume
      // `excludeSurfaceTags', which are to be removed. The remaining faces of the
      // volume become the walls of the hollowed solid, with thickness `offset'. If
      // `tag' is positive, set the tag explicitly; otherwise a new tag is selected
      // automatically.
      inline void addThickSolid(const int volumeTag,
                                const std::vector<int> & excludeSurfaceTags,
                                const double offset,
                                gmsh::vectorpair & outDimTags,
                                const int tag = -1)
      {
        int ierr = 0;
        int *api_excludeSurfaceTags_; size_t api_excludeSurfaceTags_n_; vector2ptr(excludeSurfaceTags, &api_excludeSurfaceTags_, &api_excludeSurfaceTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelOccAddThickSolid(volumeTag, api_excludeSurfaceTags_, api_excludeSurfaceTags_n_, offset, &api_outDimTags_, &api_outDimTags_n_, tag, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_excludeSurfaceTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Extrude the entities `dimTags' in the OpenCASCADE CAD representation,
      // using a translation along (`dx', `dy', `dz'). Return extruded entities in
      // `outDimTags'. If `numElements' is not empty, also extrude the mesh: the
      // entries in `numElements' give the number of elements in each layer. If
      // `height' is not empty, it provides the (cumulative) height of the
      // different layers, normalized to 1. If `recombine' is set, recombine the
      // mesh in the layers.
      inline void extrude(const gmsh::vectorpair & dimTags,
                          const double dx,
                          const double dy,
                          const double dz,
                          gmsh::vectorpair & outDimTags,
                          const std::vector<int> & numElements = std::vector<int>(),
                          const std::vector<double> & heights = std::vector<double>(),
                          const bool recombine = false)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int *api_numElements_; size_t api_numElements_n_; vector2ptr(numElements, &api_numElements_, &api_numElements_n_);
        double *api_heights_; size_t api_heights_n_; vector2ptr(heights, &api_heights_, &api_heights_n_);
        gmshModelOccExtrude(api_dimTags_, api_dimTags_n_, dx, dy, dz, &api_outDimTags_, &api_outDimTags_n_, api_numElements_, api_numElements_n_, api_heights_, api_heights_n_, (int)recombine, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        gmshFree(api_numElements_);
        gmshFree(api_heights_);
      }

      // Extrude the entities `dimTags' in the OpenCASCADE CAD representation,
      // using a rotation of `angle' radians around the axis of revolution defined
      // by the point (`x', `y', `z') and the direction (`ax', `ay', `az'). Return
      // extruded entities in `outDimTags'. If `numElements' is not empty, also
      // extrude the mesh: the entries in `numElements' give the number of elements
      // in each layer. If `height' is not empty, it provides the (cumulative)
      // height of the different layers, normalized to 1. When the mesh is extruded
      // the angle should be strictly smaller than 2*Pi. If `recombine' is set,
      // recombine the mesh in the layers.
      inline void revolve(const gmsh::vectorpair & dimTags,
                          const double x,
                          const double y,
                          const double z,
                          const double ax,
                          const double ay,
                          const double az,
                          const double angle,
                          gmsh::vectorpair & outDimTags,
                          const std::vector<int> & numElements = std::vector<int>(),
                          const std::vector<double> & heights = std::vector<double>(),
                          const bool recombine = false)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int *api_numElements_; size_t api_numElements_n_; vector2ptr(numElements, &api_numElements_, &api_numElements_n_);
        double *api_heights_; size_t api_heights_n_; vector2ptr(heights, &api_heights_, &api_heights_n_);
        gmshModelOccRevolve(api_dimTags_, api_dimTags_n_, x, y, z, ax, ay, az, angle, &api_outDimTags_, &api_outDimTags_n_, api_numElements_, api_numElements_n_, api_heights_, api_heights_n_, (int)recombine, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        gmshFree(api_numElements_);
        gmshFree(api_heights_);
      }

      // Add a pipe in the OpenCASCADE CAD representation, by extruding the
      // entities `dimTags' along the wire `wireTag'. The type of sweep can be
      // specified with `trihedron' (possible values: "DiscreteTrihedron",
      // "CorrectedFrenet", "Fixed", "Frenet", "ConstantNormal", "Darboux",
      // "GuideAC", "GuidePlan", "GuideACWithContact", "GuidePlanWithContact"). If
      // `trihedron' is not provided, "DiscreteTrihedron" is assumed. Return the
      // pipe in `outDimTags'.
      inline void addPipe(const gmsh::vectorpair & dimTags,
                          const int wireTag,
                          gmsh::vectorpair & outDimTags,
                          const std::string & trihedron = "")
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelOccAddPipe(api_dimTags_, api_dimTags_n_, wireTag, &api_outDimTags_, &api_outDimTags_n_, trihedron.c_str(), &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Fillet the volumes `volumeTags' on the curves `curveTags' with radii
      // `radii'. The `radii' vector can either contain a single radius, as many
      // radii as `curveTags', or twice as many as `curveTags' (in which case
      // different radii are provided for the begin and end points of the curves).
      // Return the filleted entities in `outDimTags'. Remove the original volume
      // if `removeVolume' is set.
      inline void fillet(const std::vector<int> & volumeTags,
                         const std::vector<int> & curveTags,
                         const std::vector<double> & radii,
                         gmsh::vectorpair & outDimTags,
                         const bool removeVolume = true)
      {
        int ierr = 0;
        int *api_volumeTags_; size_t api_volumeTags_n_; vector2ptr(volumeTags, &api_volumeTags_, &api_volumeTags_n_);
        int *api_curveTags_; size_t api_curveTags_n_; vector2ptr(curveTags, &api_curveTags_, &api_curveTags_n_);
        double *api_radii_; size_t api_radii_n_; vector2ptr(radii, &api_radii_, &api_radii_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelOccFillet(api_volumeTags_, api_volumeTags_n_, api_curveTags_, api_curveTags_n_, api_radii_, api_radii_n_, &api_outDimTags_, &api_outDimTags_n_, (int)removeVolume, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_volumeTags_);
        gmshFree(api_curveTags_);
        gmshFree(api_radii_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Chamfer the volumes `volumeTags' on the curves `curveTags' with distances
      // `distances' measured on surfaces `surfaceTags'. The `distances' vector can
      // either contain a single distance, as many distances as `curveTags' and
      // `surfaceTags', or twice as many as `curveTags' and `surfaceTags' (in which
      // case the first in each pair is measured on the corresponding surface in
      // `surfaceTags', the other on the other adjacent surface). Return the
      // chamfered entities in `outDimTags'. Remove the original volume if
      // `removeVolume' is set.
      inline void chamfer(const std::vector<int> & volumeTags,
                          const std::vector<int> & curveTags,
                          const std::vector<int> & surfaceTags,
                          const std::vector<double> & distances,
                          gmsh::vectorpair & outDimTags,
                          const bool removeVolume = true)
      {
        int ierr = 0;
        int *api_volumeTags_; size_t api_volumeTags_n_; vector2ptr(volumeTags, &api_volumeTags_, &api_volumeTags_n_);
        int *api_curveTags_; size_t api_curveTags_n_; vector2ptr(curveTags, &api_curveTags_, &api_curveTags_n_);
        int *api_surfaceTags_; size_t api_surfaceTags_n_; vector2ptr(surfaceTags, &api_surfaceTags_, &api_surfaceTags_n_);
        double *api_distances_; size_t api_distances_n_; vector2ptr(distances, &api_distances_, &api_distances_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelOccChamfer(api_volumeTags_, api_volumeTags_n_, api_curveTags_, api_curveTags_n_, api_surfaceTags_, api_surfaceTags_n_, api_distances_, api_distances_n_, &api_outDimTags_, &api_outDimTags_n_, (int)removeVolume, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_volumeTags_);
        gmshFree(api_curveTags_);
        gmshFree(api_surfaceTags_);
        gmshFree(api_distances_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Compute the boolean union (the fusion) of the entities `objectDimTags' and
      // `toolDimTags' in the OpenCASCADE CAD representation. Return the resulting
      // entities in `outDimTags'. If `tag' is positive, try to set the tag
      // explicitly (only valid if the boolean operation results in a single
      // entity). Remove the object if `removeObject' is set. Remove the tool if
      // `removeTool' is set.
      inline void fuse(const gmsh::vectorpair & objectDimTags,
                       const gmsh::vectorpair & toolDimTags,
                       gmsh::vectorpair & outDimTags,
                       std::vector<gmsh::vectorpair> & outDimTagsMap,
                       const int tag = -1,
                       const bool removeObject = true,
                       const bool removeTool = true)
      {
        int ierr = 0;
        int *api_objectDimTags_; size_t api_objectDimTags_n_; vectorpair2intptr(objectDimTags, &api_objectDimTags_, &api_objectDimTags_n_);
        int *api_toolDimTags_; size_t api_toolDimTags_n_; vectorpair2intptr(toolDimTags, &api_toolDimTags_, &api_toolDimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int **api_outDimTagsMap_; size_t *api_outDimTagsMap_n_, api_outDimTagsMap_nn_;
        gmshModelOccFuse(api_objectDimTags_, api_objectDimTags_n_, api_toolDimTags_, api_toolDimTags_n_, &api_outDimTags_, &api_outDimTags_n_, &api_outDimTagsMap_, &api_outDimTagsMap_n_, &api_outDimTagsMap_nn_, tag, (int)removeObject, (int)removeTool, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_objectDimTags_);
        gmshFree(api_toolDimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        outDimTagsMap.resize(api_outDimTagsMap_nn_); for(size_t i = 0; i < api_outDimTagsMap_nn_; ++i){ outDimTagsMap[i].resize(api_outDimTagsMap_n_[i] / 2); for(size_t j = 0; j < api_outDimTagsMap_n_[i] / 2; ++j){ outDimTagsMap[i][j].first = api_outDimTagsMap_[i][j * 2 + 0]; outDimTagsMap[i][j].second = api_outDimTagsMap_[i][j * 2 + 1]; } gmshFree(api_outDimTagsMap_[i]); } gmshFree(api_outDimTagsMap_); gmshFree(api_outDimTagsMap_n_);
      }

      // Compute the boolean intersection (the common parts) of the entities
      // `objectDimTags' and `toolDimTags' in the OpenCASCADE CAD representation.
      // Return the resulting entities in `outDimTags'. If `tag' is positive, try
      // to set the tag explicitly (only valid if the boolean operation results in
      // a single entity). Remove the object if `removeObject' is set. Remove the
      // tool if `removeTool' is set.
      inline void intersect(const gmsh::vectorpair & objectDimTags,
                            const gmsh::vectorpair & toolDimTags,
                            gmsh::vectorpair & outDimTags,
                            std::vector<gmsh::vectorpair> & outDimTagsMap,
                            const int tag = -1,
                            const bool removeObject = true,
                            const bool removeTool = true)
      {
        int ierr = 0;
        int *api_objectDimTags_; size_t api_objectDimTags_n_; vectorpair2intptr(objectDimTags, &api_objectDimTags_, &api_objectDimTags_n_);
        int *api_toolDimTags_; size_t api_toolDimTags_n_; vectorpair2intptr(toolDimTags, &api_toolDimTags_, &api_toolDimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int **api_outDimTagsMap_; size_t *api_outDimTagsMap_n_, api_outDimTagsMap_nn_;
        gmshModelOccIntersect(api_objectDimTags_, api_objectDimTags_n_, api_toolDimTags_, api_toolDimTags_n_, &api_outDimTags_, &api_outDimTags_n_, &api_outDimTagsMap_, &api_outDimTagsMap_n_, &api_outDimTagsMap_nn_, tag, (int)removeObject, (int)removeTool, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_objectDimTags_);
        gmshFree(api_toolDimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        outDimTagsMap.resize(api_outDimTagsMap_nn_); for(size_t i = 0; i < api_outDimTagsMap_nn_; ++i){ outDimTagsMap[i].resize(api_outDimTagsMap_n_[i] / 2); for(size_t j = 0; j < api_outDimTagsMap_n_[i] / 2; ++j){ outDimTagsMap[i][j].first = api_outDimTagsMap_[i][j * 2 + 0]; outDimTagsMap[i][j].second = api_outDimTagsMap_[i][j * 2 + 1]; } gmshFree(api_outDimTagsMap_[i]); } gmshFree(api_outDimTagsMap_); gmshFree(api_outDimTagsMap_n_);
      }

      // Compute the boolean difference between the entities `objectDimTags' and
      // `toolDimTags' in the OpenCASCADE CAD representation. Return the resulting
      // entities in `outDimTags'. If `tag' is positive, try to set the tag
      // explicitly (only valid if the boolean operation results in a single
      // entity). Remove the object if `removeObject' is set. Remove the tool if
      // `removeTool' is set.
      inline void cut(const gmsh::vectorpair & objectDimTags,
                      const gmsh::vectorpair & toolDimTags,
                      gmsh::vectorpair & outDimTags,
                      std::vector<gmsh::vectorpair> & outDimTagsMap,
                      const int tag = -1,
                      const bool removeObject = true,
                      const bool removeTool = true)
      {
        int ierr = 0;
        int *api_objectDimTags_; size_t api_objectDimTags_n_; vectorpair2intptr(objectDimTags, &api_objectDimTags_, &api_objectDimTags_n_);
        int *api_toolDimTags_; size_t api_toolDimTags_n_; vectorpair2intptr(toolDimTags, &api_toolDimTags_, &api_toolDimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int **api_outDimTagsMap_; size_t *api_outDimTagsMap_n_, api_outDimTagsMap_nn_;
        gmshModelOccCut(api_objectDimTags_, api_objectDimTags_n_, api_toolDimTags_, api_toolDimTags_n_, &api_outDimTags_, &api_outDimTags_n_, &api_outDimTagsMap_, &api_outDimTagsMap_n_, &api_outDimTagsMap_nn_, tag, (int)removeObject, (int)removeTool, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_objectDimTags_);
        gmshFree(api_toolDimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        outDimTagsMap.resize(api_outDimTagsMap_nn_); for(size_t i = 0; i < api_outDimTagsMap_nn_; ++i){ outDimTagsMap[i].resize(api_outDimTagsMap_n_[i] / 2); for(size_t j = 0; j < api_outDimTagsMap_n_[i] / 2; ++j){ outDimTagsMap[i][j].first = api_outDimTagsMap_[i][j * 2 + 0]; outDimTagsMap[i][j].second = api_outDimTagsMap_[i][j * 2 + 1]; } gmshFree(api_outDimTagsMap_[i]); } gmshFree(api_outDimTagsMap_); gmshFree(api_outDimTagsMap_n_);
      }

      // Compute the boolean fragments (general fuse) resulting from the
      // intersection of the entities `objectDimTags' and `toolDimTags' in the
      // OpenCASCADE CAD representation, making all iterfaces conformal. When
      // applied to entities of different dimensions, the lower dimensional
      // entities will be automatically embedded in the higher dimensional entities
      // if they are not on their boundary. Return the resulting entities in
      // `outDimTags'. If `tag' is positive, try to set the tag explicitly (only
      // valid if the boolean operation results in a single entity). Remove the
      // object if `removeObject' is set. Remove the tool if `removeTool' is set.
      inline void fragment(const gmsh::vectorpair & objectDimTags,
                           const gmsh::vectorpair & toolDimTags,
                           gmsh::vectorpair & outDimTags,
                           std::vector<gmsh::vectorpair> & outDimTagsMap,
                           const int tag = -1,
                           const bool removeObject = true,
                           const bool removeTool = true)
      {
        int ierr = 0;
        int *api_objectDimTags_; size_t api_objectDimTags_n_; vectorpair2intptr(objectDimTags, &api_objectDimTags_, &api_objectDimTags_n_);
        int *api_toolDimTags_; size_t api_toolDimTags_n_; vectorpair2intptr(toolDimTags, &api_toolDimTags_, &api_toolDimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int **api_outDimTagsMap_; size_t *api_outDimTagsMap_n_, api_outDimTagsMap_nn_;
        gmshModelOccFragment(api_objectDimTags_, api_objectDimTags_n_, api_toolDimTags_, api_toolDimTags_n_, &api_outDimTags_, &api_outDimTags_n_, &api_outDimTagsMap_, &api_outDimTagsMap_n_, &api_outDimTagsMap_nn_, tag, (int)removeObject, (int)removeTool, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_objectDimTags_);
        gmshFree(api_toolDimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        outDimTagsMap.resize(api_outDimTagsMap_nn_); for(size_t i = 0; i < api_outDimTagsMap_nn_; ++i){ outDimTagsMap[i].resize(api_outDimTagsMap_n_[i] / 2); for(size_t j = 0; j < api_outDimTagsMap_n_[i] / 2; ++j){ outDimTagsMap[i][j].first = api_outDimTagsMap_[i][j * 2 + 0]; outDimTagsMap[i][j].second = api_outDimTagsMap_[i][j * 2 + 1]; } gmshFree(api_outDimTagsMap_[i]); } gmshFree(api_outDimTagsMap_); gmshFree(api_outDimTagsMap_n_);
      }

      // Translate the entities `dimTags' in the OpenCASCADE CAD representation
      // along (`dx', `dy', `dz').
      inline void translate(const gmsh::vectorpair & dimTags,
                            const double dx,
                            const double dy,
                            const double dz)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelOccTranslate(api_dimTags_, api_dimTags_n_, dx, dy, dz, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Rotate the entities `dimTags' in the OpenCASCADE CAD representation by
      // `angle' radians around the axis of revolution defined by the point (`x',
      // `y', `z') and the direction (`ax', `ay', `az').
      inline void rotate(const gmsh::vectorpair & dimTags,
                         const double x,
                         const double y,
                         const double z,
                         const double ax,
                         const double ay,
                         const double az,
                         const double angle)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelOccRotate(api_dimTags_, api_dimTags_n_, x, y, z, ax, ay, az, angle, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Scale the entities `dimTags' in the OpenCASCADE CAD representation by
      // factors `a', `b' and `c' along the three coordinate axes; use (`x', `y',
      // `z') as the center of the homothetic transformation.
      inline void dilate(const gmsh::vectorpair & dimTags,
                         const double x,
                         const double y,
                         const double z,
                         const double a,
                         const double b,
                         const double c)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelOccDilate(api_dimTags_, api_dimTags_n_, x, y, z, a, b, c, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Mirror the entities `dimTags' in the OpenCASCADE CAD representation, with
      // respect to the plane of equation `a' * x + `b' * y + `c' * z + `d' = 0.
      inline void mirror(const gmsh::vectorpair & dimTags,
                         const double a,
                         const double b,
                         const double c,
                         const double d)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelOccMirror(api_dimTags_, api_dimTags_n_, a, b, c, d, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Mirror the entities `dimTags' in the OpenCASCADE CAD representation, with
      // respect to the plane of equation `a' * x + `b' * y + `c' * z + `d' = 0.
      // (This is a deprecated synonym for `mirror'.)
      inline void symmetrize(const gmsh::vectorpair & dimTags,
                             const double a,
                             const double b,
                             const double c,
                             const double d)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelOccSymmetrize(api_dimTags_, api_dimTags_n_, a, b, c, d, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Apply a general affine transformation matrix `affineTransform' (16 entries
      // of a 4x4 matrix, by row; only the 12 first can be provided for
      // convenience) to the entities `dimTags' in the OpenCASCADE CAD
      // representation.
      inline void affineTransform(const gmsh::vectorpair & dimTags,
                                  const std::vector<double> & affineTransform)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        double *api_affineTransform_; size_t api_affineTransform_n_; vector2ptr(affineTransform, &api_affineTransform_, &api_affineTransform_n_);
        gmshModelOccAffineTransform(api_dimTags_, api_dimTags_n_, api_affineTransform_, api_affineTransform_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        gmshFree(api_affineTransform_);
      }

      // Copy the entities `dimTags' in the OpenCASCADE CAD representation; the new
      // entities are returned in `outDimTags'.
      inline void copy(const gmsh::vectorpair & dimTags,
                       gmsh::vectorpair & outDimTags)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelOccCopy(api_dimTags_, api_dimTags_n_, &api_outDimTags_, &api_outDimTags_n_, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Remove the entities `dimTags' in the OpenCASCADE CAD representation,
      // provided that they are not on the boundary of higher-dimensional entities.
      // If `recursive' is true, remove all the entities on their boundaries, down
      // to dimension 0.
      inline void remove(const gmsh::vectorpair & dimTags,
                         const bool recursive = false)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelOccRemove(api_dimTags_, api_dimTags_n_, (int)recursive, &ierr);
        if(ierr) throwLastError();
        gmshFree(api_dimTags_);
      }

      // Remove all duplicate entities in the OpenCASCADE CAD representation
      // (different entities at the same geometrical location) after intersecting
      // (using boolean fragments) all highest dimensional entities.
      inline void removeAllDuplicates()
      {
        int ierr = 0;
        gmshModelOccRemoveAllDuplicates(&ierr);
        if(ierr) throwLastError();
      }

      // Apply various healing procedures to the entities `dimTags' (or to all the
      // entities in the model if `dimTags' is empty) in the OpenCASCADE CAD
      // representation. Return the healed entities in `outDimTags'. Available
      // healing options are listed in the Gmsh reference manual.
      inline void healShapes(gmsh::vectorpair & outDimTags,
                             const gmsh::vectorpair & dimTags = gmsh::vectorpair(),
                             const double tolerance = 1e-8,
                             const bool fixDegenerated = true,
                             const bool fixSmallEdges = true,
                             const bool fixSmallFaces = true,
                             const bool sewFaces = true,
                             const bool makeSolids = true)
      {
        int ierr = 0;
        int *api_outDimTags_; size_t api_outDimTags_n_;
        int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
        gmshModelOccHealShapes(&api_outDimTags_, &api_outDimTags_n_, api_dimTags_, api_dimTags_n_, tolerance, (int)fixDegenerated, (int)fixSmallEdges, (int)fixSmallFaces, (int)sewFaces, (int)makeSolids, &ierr);
        if(ierr) throwLastError();
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
        gmshFree(api_dimTags_);
      }

      // Import BREP, STEP or IGES shapes from the file `fileName' in the
      // OpenCASCADE CAD representation. The imported entities are returned in
      // `outDimTags'. If the optional argument `highestDimOnly' is set, only
      // import the highest dimensional entities in the file. The optional argument
      // `format' can be used to force the format of the file (currently "brep",
      // "step" or "iges").
      inline void importShapes(const std::string & fileName,
                               gmsh::vectorpair & outDimTags,
                               const bool highestDimOnly = true,
                               const std::string & format = "")
      {
        int ierr = 0;
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelOccImportShapes(fileName.c_str(), &api_outDimTags_, &api_outDimTags_n_, (int)highestDimOnly, format.c_str(), &ierr);
        if(ierr) throwLastError();
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Imports an OpenCASCADE `shape' by providing a pointer to a native
      // OpenCASCADE `TopoDS_Shape' object (passed as a pointer to void). The
      // imported entities are returned in `outDimTags'. If the optional argument
      // `highestDimOnly' is set, only import the highest dimensional entities in
      // `shape'. In Python, this function can be used for integration with
      // PythonOCC, in which the SwigPyObject pointer of `TopoDS_Shape' must be
      // passed as an int to `shape', i.e., `shape = int(pythonocc_shape.this)'.
      // Warning: this function is unsafe, as providing an invalid pointer will
      // lead to undefined behavior.
      inline void importShapesNativePointer(const void * shape,
                                            gmsh::vectorpair & outDimTags,
                                            const bool highestDimOnly = true)
      {
        int ierr = 0;
        int *api_outDimTags_; size_t api_outDimTags_n_;
        gmshModelOccImportShapesNativePointer(shape, &api_outDimTags_, &api_outDimTags_n_, (int)highestDimOnly, &ierr);
        if(ierr) throwLastError();
        outDimTags.resize(api_outDimTags_n_ / 2); for(size_t i = 0; i < api_outDimTags_n_ / 2; ++i){ outDimTags[i].first = api_outDimTags_[i * 2 + 0]; outDimTags[i].second = api_outDimTags_[i * 2 + 1]; } gmshFree(api_outDimTags_);
      }

      // Get all the OpenCASCADE entities. If `dim' is >= 0, return only the
      // entities of the specified dimension (e.g. points if `dim' == 0). The
      // entities are returned as a vector of (dim, tag) integer pairs.
      inline void getEntities(gmsh::vectorpair & dimTags,
                              const int dim = -1)
      {
        int ierr = 0;
        int *api_dimTags_; size_t api_dimTags_n_;
        gmshModelOccGetEntities(&api_dimTags_, &api_dimTags_n_, dim, &ierr);
        if(ierr) throwLastError();
        dimTags.resize(api_dimTags_n_ / 2); for(size_t i = 0; i < api_dimTags_n_ / 2; ++i){ dimTags[i].first = api_dimTags_[i * 2 + 0]; dimTags[i].second = api_dimTags_[i * 2 + 1]; } gmshFree(api_dimTags_);
      }

      // Get the OpenCASCADE entities in the bounding box defined by the two points
      // (`xmin', `ymin', `zmin') and (`xmax', `ymax', `zmax'). If `dim' is >= 0,
      // return only the entities of the specified dimension (e.g. points if `dim'
      // == 0).
      inline void getEntitiesInBoundingBox(const double xmin,
                                           const double ymin,
                                           const double zmin,
                                           const double xmax,
                                           const double ymax,
                                           const double zmax,
                                           gmsh::vectorpair & tags,
                                           const int dim = -1)
      {
        int ierr = 0;
        int *api_tags_; size_t api_tags_n_;
        gmshModelOccGetEntitiesInBoundingBox(xmin, ymin, zmin, xmax, ymax, zmax, &api_tags_, &api_tags_n_, dim, &ierr);
        if(ierr) throwLastError();
        tags.resize(api_tags_n_ / 2); for(size_t i = 0; i < api_tags_n_ / 2; ++i){ tags[i].first = api_tags_[i * 2 + 0]; tags[i].second = api_tags_[i * 2 + 1]; } gmshFree(api_tags_);
      }

      // Get the bounding box (`xmin', `ymin', `zmin'), (`xmax', `ymax', `zmax') of
      // the OpenCASCADE entity of dimension `dim' and tag `tag'.
      inline void getBoundingBox(const int dim,
                                 const int tag,
                                 double & xmin,
                                 double & ymin,
                                 double & zmin,
                                 double & xmax,
                                 double & ymax,
                                 double & zmax)
      {
        int ierr = 0;
        gmshModelOccGetBoundingBox(dim, tag, &xmin, &ymin, &zmin, &xmax, &ymax, &zmax, &ierr);
        if(ierr) throwLastError();
      }

      // Get the mass of the OpenCASCADE entity of dimension `dim' and tag `tag'.
      inline void getMass(const int dim,
                          const int tag,
                          double & mass)
      {
        int ierr = 0;
        gmshModelOccGetMass(dim, tag, &mass, &ierr);
        if(ierr) throwLastError();
      }

      // Get the center of mass of the OpenCASCADE entity of dimension `dim' and
      // tag `tag'.
      inline void getCenterOfMass(const int dim,
                                  const int tag,
                                  double & x,
                                  double & y,
                                  double & z)
      {
        int ierr = 0;
        gmshModelOccGetCenterOfMass(dim, tag, &x, &y, &z, &ierr);
        if(ierr) throwLastError();
      }

      // Get the matrix of inertia (by row) of the OpenCASCADE entity of dimension
      // `dim' and tag `tag'.
      inline void getMatrixOfInertia(const int dim,
                                     const int tag,
                                     std::vector<double> & mat)
      {
        int ierr = 0;
        double *api_mat_; size_t api_mat_n_;
        gmshModelOccGetMatrixOfInertia(dim, tag, &api_mat_, &api_mat_n_, &ierr);
        if(ierr) throwLastError();
        mat.assign(api_mat_, api_mat_ + api_mat_n_); gmshFree(api_mat_);
      }

      // Get the maximum tag of entities of dimension `dim' in the OpenCASCADE CAD
      // representation.
      inline int getMaxTag(const int dim)
      {
        int ierr = 0;
        int result_api_ = gmshModelOccGetMaxTag(dim, &ierr);
        if(ierr) throwLastError();
        return result_api_;
      }

      // Set the maximum tag `maxTag' for entities of dimension `dim' in the
      // OpenCASCADE CAD representation.
      inline void setMaxTag(const int dim,
                            const int maxTag)
      {
        int ierr = 0;
        gmshModelOccSetMaxTag(dim, maxTag, &ierr);
        if(ierr) throwLastError();
      }

      // Synchronize the OpenCASCADE CAD representation with the current Gmsh
      // model. This can be called at any time, but since it involves a non trivial
      // amount of processing, the number of synchronization points should normally
      // be minimized. Without synchronization the entities in the OpenCASCADE CAD
      // representation are not available to any function outside of the
      // OpenCASCADE CAD kernel functions.
      inline void synchronize()
      {
        int ierr = 0;
        gmshModelOccSynchronize(&ierr);
        if(ierr) throwLastError();
      }

      namespace mesh { // OpenCASCADE CAD kernel meshing constraints

        // Set a mesh size constraint on the entities `dimTags' in the OpenCASCADE
        // CAD representation. Currently only entities of dimension 0 (points) are
        // handled.
        inline void setSize(const gmsh::vectorpair & dimTags,
                            const double size)
        {
          int ierr = 0;
          int *api_dimTags_; size_t api_dimTags_n_; vectorpair2intptr(dimTags, &api_dimTags_, &api_dimTags_n_);
          gmshModelOccMeshSetSize(api_dimTags_, api_dimTags_n_, size, &ierr);
          if(ierr) throwLastError();
          gmshFree(api_dimTags_);
        }

      } // namespace mesh

    } // namespace occ

  } // namespace model

  namespace view { // Post-processing view functions

    // Add a new post-processing view, with name `name'. If `tag' is positive use
    // it (and remove the view with that tag if it already exists), otherwise
    // associate a new tag. Return the view tag.
    inline int add(const std::string & name,
                   const int tag = -1)
    {
      int ierr = 0;
      int result_api_ = gmshViewAdd(name.c_str(), tag, &ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

    // Remove the view with tag `tag'.
    inline void remove(const int tag)
    {
      int ierr = 0;
      gmshViewRemove(tag, &ierr);
      if(ierr) throwLastError();
    }

    // Get the index of the view with tag `tag' in the list of currently loaded
    // views. This dynamic index (it can change when views are removed) is used to
    // access view options.
    inline int getIndex(const int tag)
    {
      int ierr = 0;
      int result_api_ = gmshViewGetIndex(tag, &ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

    // Get the tags of all views.
    inline void getTags(std::vector<int> & tags)
    {
      int ierr = 0;
      int *api_tags_; size_t api_tags_n_;
      gmshViewGetTags(&api_tags_, &api_tags_n_, &ierr);
      if(ierr) throwLastError();
      tags.assign(api_tags_, api_tags_ + api_tags_n_); gmshFree(api_tags_);
    }

    // Add model-based post-processing data to the view with tag `tag'. `modelName'
    // identifies the model the data is attached to. `dataType' specifies the type
    // of data, currently either "NodeData", "ElementData" or "ElementNodeData".
    // `step' specifies the identifier (>= 0) of the data in a sequence. `tags'
    // gives the tags of the nodes or elements in the mesh to which the data is
    // associated. `data' is a vector of the same length as `tags': each entry is
    // the vector of double precision numbers representing the data associated with
    // the corresponding tag. The optional `time' argument associate a time value
    // with the data. `numComponents' gives the number of data components (1 for
    // scalar data, 3 for vector data, etc.) per entity; if negative, it is
    // automatically inferred (when possible) from the input data. `partition'
    // allows one to specify data in several sub-sets.
    inline void addModelData(const int tag,
                             const int step,
                             const std::string & modelName,
                             const std::string & dataType,
                             const std::vector<std::size_t> & tags,
                             const std::vector<std::vector<double> > & data,
                             const double time = 0.,
                             const int numComponents = -1,
                             const int partition = 0)
    {
      int ierr = 0;
      size_t *api_tags_; size_t api_tags_n_; vector2ptr(tags, &api_tags_, &api_tags_n_);
      double **api_data_; size_t *api_data_n_, api_data_nn_; vectorvector2ptrptr(data, &api_data_, &api_data_n_, &api_data_nn_);
      gmshViewAddModelData(tag, step, modelName.c_str(), dataType.c_str(), api_tags_, api_tags_n_, (const double **)api_data_, api_data_n_, api_data_nn_, time, numComponents, partition, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_tags_);
      for(size_t i = 0; i < api_data_nn_; ++i){ gmshFree(api_data_[i]); } gmshFree(api_data_); gmshFree(api_data_n_);
    }

    // Add homogeneous model-based post-processing data to the view with tag `tag'.
    // The arguments have the same meaning as in `addModelData', except that `data'
    // is supposed to be homogeneous and is thus flattened in a single vector. For
    // data types that can lead to different data sizes per tag (like
    // "ElementNodeData"), the data should be padded.
    inline void addHomogeneousModelData(const int tag,
                                        const int step,
                                        const std::string & modelName,
                                        const std::string & dataType,
                                        const std::vector<std::size_t> & tags,
                                        const std::vector<double> & data,
                                        const double time = 0.,
                                        const int numComponents = -1,
                                        const int partition = 0)
    {
      int ierr = 0;
      size_t *api_tags_; size_t api_tags_n_; vector2ptr(tags, &api_tags_, &api_tags_n_);
      double *api_data_; size_t api_data_n_; vector2ptr(data, &api_data_, &api_data_n_);
      gmshViewAddHomogeneousModelData(tag, step, modelName.c_str(), dataType.c_str(), api_tags_, api_tags_n_, api_data_, api_data_n_, time, numComponents, partition, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_tags_);
      gmshFree(api_data_);
    }

    // Get model-based post-processing data from the view with tag `tag' at step
    // `step'. Return the `data' associated to the nodes or the elements with tags
    // `tags', as well as the `dataType' and the number of components
    // `numComponents'.
    inline void getModelData(const int tag,
                             const int step,
                             std::string & dataType,
                             std::vector<std::size_t> & tags,
                             std::vector<std::vector<double> > & data,
                             double & time,
                             int & numComponents)
    {
      int ierr = 0;
      char *api_dataType_;
      size_t *api_tags_; size_t api_tags_n_;
      double **api_data_; size_t *api_data_n_, api_data_nn_;
      gmshViewGetModelData(tag, step, &api_dataType_, &api_tags_, &api_tags_n_, &api_data_, &api_data_n_, &api_data_nn_, &time, &numComponents, &ierr);
      if(ierr) throwLastError();
      dataType = std::string(api_dataType_); gmshFree(api_dataType_);
      tags.assign(api_tags_, api_tags_ + api_tags_n_); gmshFree(api_tags_);
      data.resize(api_data_nn_); for(size_t i = 0; i < api_data_nn_; ++i){ data[i].assign(api_data_[i], api_data_[i] + api_data_n_[i]); gmshFree(api_data_[i]); } gmshFree(api_data_); gmshFree(api_data_n_);
    }

    // Get homogeneous model-based post-processing data from the view with tag
    // `tag' at step `step'. The arguments have the same meaning as in
    // `getModelData', except that `data' is returned flattened in a single vector,
    // with the appropriate padding if necessary.
    inline void getHomogeneousModelData(const int tag,
                                        const int step,
                                        std::string & dataType,
                                        std::vector<std::size_t> & tags,
                                        std::vector<double> & data,
                                        double & time,
                                        int & numComponents)
    {
      int ierr = 0;
      char *api_dataType_;
      size_t *api_tags_; size_t api_tags_n_;
      double *api_data_; size_t api_data_n_;
      gmshViewGetHomogeneousModelData(tag, step, &api_dataType_, &api_tags_, &api_tags_n_, &api_data_, &api_data_n_, &time, &numComponents, &ierr);
      if(ierr) throwLastError();
      dataType = std::string(api_dataType_); gmshFree(api_dataType_);
      tags.assign(api_tags_, api_tags_ + api_tags_n_); gmshFree(api_tags_);
      data.assign(api_data_, api_data_ + api_data_n_); gmshFree(api_data_);
    }

    // Add list-based post-processing data to the view with tag `tag'. List-based
    // datasets are independent from any model and any mesh. `dataType' identifies
    // the data by concatenating the field type ("S" for scalar, "V" for vector,
    // "T" for tensor) and the element type ("P" for point, "L" for line, "T" for
    // triangle, "S" for tetrahedron, "I" for prism, "H" for hexaHedron, "Y" for
    // pyramid). For example `dataType' should be "ST" for a scalar field on
    // triangles. `numEle' gives the number of elements in the data. `data'
    // contains the data for the `numEle' elements, concatenated, with node
    // coordinates followed by values per node, repeated for each step: [e1x1, ...,
    // e1xn, e1y1, ..., e1yn, e1z1, ..., e1zn, e1v1..., e1vN, e2x1, ...].
    inline void addListData(const int tag,
                            const std::string & dataType,
                            const int numEle,
                            const std::vector<double> & data)
    {
      int ierr = 0;
      double *api_data_; size_t api_data_n_; vector2ptr(data, &api_data_, &api_data_n_);
      gmshViewAddListData(tag, dataType.c_str(), numEle, api_data_, api_data_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_data_);
    }

    // Get list-based post-processing data from the view with tag `tag'. Return the
    // types `dataTypes', the number of elements `numElements' for each data type
    // and the `data' for each data type.
    inline void getListData(const int tag,
                            std::vector<std::string> & dataType,
                            std::vector<int> & numElements,
                            std::vector<std::vector<double> > & data)
    {
      int ierr = 0;
      char **api_dataType_; size_t api_dataType_n_;
      int *api_numElements_; size_t api_numElements_n_;
      double **api_data_; size_t *api_data_n_, api_data_nn_;
      gmshViewGetListData(tag, &api_dataType_, &api_dataType_n_, &api_numElements_, &api_numElements_n_, &api_data_, &api_data_n_, &api_data_nn_, &ierr);
      if(ierr) throwLastError();
      dataType.resize(api_dataType_n_); for(size_t i = 0; i < api_dataType_n_; ++i){ dataType[i] = std::string(api_dataType_[i]); gmshFree(api_dataType_[i]); } gmshFree(api_dataType_);
      numElements.assign(api_numElements_, api_numElements_ + api_numElements_n_); gmshFree(api_numElements_);
      data.resize(api_data_nn_); for(size_t i = 0; i < api_data_nn_; ++i){ data[i].assign(api_data_[i], api_data_[i] + api_data_n_[i]); gmshFree(api_data_[i]); } gmshFree(api_data_); gmshFree(api_data_n_);
    }

    // Add a string to a list-based post-processing view with tag `tag'. If `coord'
    // contains 3 coordinates the string is positioned in the 3D model space ("3D
    // string"); if it contains 2 coordinates it is positioned in the 2D graphics
    // viewport ("2D string"). `data' contains one or more (for multistep views)
    // strings. `style' contains key-value pairs of styling parameters,
    // concatenated. Available keys are "Font" (possible values: "Times-Roman",
    // "Times-Bold", "Times-Italic", "Times-BoldItalic", "Helvetica", "Helvetica-
    // Bold", "Helvetica-Oblique", "Helvetica-BoldOblique", "Courier", "Courier-
    // Bold", "Courier-Oblique", "Courier-BoldOblique", "Symbol", "ZapfDingbats",
    // "Screen"), "FontSize" and "Align" (possible values: "Left" or "BottomLeft",
    // "Center" or "BottomCenter", "Right" or "BottomRight", "TopLeft",
    // "TopCenter", "TopRight", "CenterLeft", "CenterCenter", "CenterRight").
    inline void addListDataString(const int tag,
                                  const std::vector<double> & coord,
                                  const std::vector<std::string> & data,
                                  const std::vector<std::string> & style = std::vector<std::string>())
    {
      int ierr = 0;
      double *api_coord_; size_t api_coord_n_; vector2ptr(coord, &api_coord_, &api_coord_n_);
      char **api_data_; size_t api_data_n_; vectorstring2charptrptr(data, &api_data_, &api_data_n_);
      char **api_style_; size_t api_style_n_; vectorstring2charptrptr(style, &api_style_, &api_style_n_);
      gmshViewAddListDataString(tag, api_coord_, api_coord_n_, api_data_, api_data_n_, api_style_, api_style_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_coord_);
      for(size_t i = 0; i < api_data_n_; ++i){ gmshFree(api_data_[i]); } gmshFree(api_data_);
      for(size_t i = 0; i < api_style_n_; ++i){ gmshFree(api_style_[i]); } gmshFree(api_style_);
    }

    // Get list-based post-processing data strings (2D strings if `dim' = 2, 3D
    // strings if `dim' = 3) from the view with tag `tag'. Return the coordinates
    // in `coord', the strings in `data' and the styles in `style'.
    inline void getListDataStrings(const int tag,
                                   const int dim,
                                   std::vector<double> & coord,
                                   std::vector<std::string> & data,
                                   std::vector<std::string> & style)
    {
      int ierr = 0;
      double *api_coord_; size_t api_coord_n_;
      char **api_data_; size_t api_data_n_;
      char **api_style_; size_t api_style_n_;
      gmshViewGetListDataStrings(tag, dim, &api_coord_, &api_coord_n_, &api_data_, &api_data_n_, &api_style_, &api_style_n_, &ierr);
      if(ierr) throwLastError();
      coord.assign(api_coord_, api_coord_ + api_coord_n_); gmshFree(api_coord_);
      data.resize(api_data_n_); for(size_t i = 0; i < api_data_n_; ++i){ data[i] = std::string(api_data_[i]); gmshFree(api_data_[i]); } gmshFree(api_data_);
      style.resize(api_style_n_); for(size_t i = 0; i < api_style_n_; ++i){ style[i] = std::string(api_style_[i]); gmshFree(api_style_[i]); } gmshFree(api_style_);
    }

    // Set interpolation matrices for the element family `type' ("Line",
    // "Triangle", "Quadrangle", "Tetrahedron", "Hexahedron", "Prism", "Pyramid")
    // in the view `tag'. The approximation of the values over an element is
    // written as a linear combination of `d' basis functions f_i(u, v, w) = sum_(j
    // = 0, ..., `d' - 1) `coef'[i][j] u^`exp'[j][0] v^`exp'[j][1] w^`exp'[j][2], i
    // = 0, ..., `d'-1, with u, v, w the coordinates in the reference element. The
    // `coef' matrix (of size `d' x `d') and the `exp' matrix (of size `d' x 3) are
    // stored as vectors, by row. If `dGeo' is positive, use `coefGeo' and `expGeo'
    // to define the interpolation of the x, y, z coordinates of the element in
    // terms of the u, v, w coordinates, in exactly the same way. If `d' < 0,
    // remove the interpolation matrices.
    inline void setInterpolationMatrices(const int tag,
                                         const std::string & type,
                                         const int d,
                                         const std::vector<double> & coef,
                                         const std::vector<double> & exp,
                                         const int dGeo = 0,
                                         const std::vector<double> & coefGeo = std::vector<double>(),
                                         const std::vector<double> & expGeo = std::vector<double>())
    {
      int ierr = 0;
      double *api_coef_; size_t api_coef_n_; vector2ptr(coef, &api_coef_, &api_coef_n_);
      double *api_exp_; size_t api_exp_n_; vector2ptr(exp, &api_exp_, &api_exp_n_);
      double *api_coefGeo_; size_t api_coefGeo_n_; vector2ptr(coefGeo, &api_coefGeo_, &api_coefGeo_n_);
      double *api_expGeo_; size_t api_expGeo_n_; vector2ptr(expGeo, &api_expGeo_, &api_expGeo_n_);
      gmshViewSetInterpolationMatrices(tag, type.c_str(), d, api_coef_, api_coef_n_, api_exp_, api_exp_n_, dGeo, api_coefGeo_, api_coefGeo_n_, api_expGeo_, api_expGeo_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_coef_);
      gmshFree(api_exp_);
      gmshFree(api_coefGeo_);
      gmshFree(api_expGeo_);
    }

    // Add a post-processing view as an `alias' of the reference view with tag
    // `refTag'. If `copyOptions' is set, copy the options of the reference view.
    // If `tag' is positive use it (and remove the view with that tag if it already
    // exists), otherwise associate a new tag. Return the view tag.
    inline int addAlias(const int refTag,
                        const bool copyOptions = false,
                        const int tag = -1)
    {
      int ierr = 0;
      int result_api_ = gmshViewAddAlias(refTag, (int)copyOptions, tag, &ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

    // Combine elements (if `what' == "elements") or steps (if `what' == "steps")
    // of all views (`how' == "all"), all visible views (`how' == "visible") or all
    // views having the same name (`how' == "name"). Remove original views if
    // `remove' is set.
    inline void combine(const std::string & what,
                        const std::string & how,
                        const bool remove = true,
                        const bool copyOptions = true)
    {
      int ierr = 0;
      gmshViewCombine(what.c_str(), how.c_str(), (int)remove, (int)copyOptions, &ierr);
      if(ierr) throwLastError();
    }

    // Probe the view `tag' for its `value' at point (`x', `y', `z'). If no match
    // is found, `value' is returned empty. Return only the value at step `step' is
    // `step' is positive. Return only values with `numComp' if `numComp' is
    // positive. Return the gradient of the `value' if `gradient' is set. If
    // `distanceMax' is zero, only return a result if an exact match inside an
    // element in the view is found; if `distanceMax' is positive and an exact
    // match is not found, return the value at the closest node if it is closer
    // than `distanceMax'; if `distanceMax' is negative and an exact match is not
    // found, always return the value at the closest node. The distance to the
    // match is returned in `distance'. Return the result from the element
    // described by its coordinates if `xElementCoord', `yElementCoord' and
    // `zElementCoord' are provided. If `dim' is >= 0, return only matches from
    // elements of the specified dimension.
    inline void probe(const int tag,
                      const double x,
                      const double y,
                      const double z,
                      std::vector<double> & value,
                      double & distance,
                      const int step = -1,
                      const int numComp = -1,
                      const bool gradient = false,
                      const double distanceMax = 0.,
                      const std::vector<double> & xElemCoord = std::vector<double>(),
                      const std::vector<double> & yElemCoord = std::vector<double>(),
                      const std::vector<double> & zElemCoord = std::vector<double>(),
                      const int dim = -1)
    {
      int ierr = 0;
      double *api_value_; size_t api_value_n_;
      double *api_xElemCoord_; size_t api_xElemCoord_n_; vector2ptr(xElemCoord, &api_xElemCoord_, &api_xElemCoord_n_);
      double *api_yElemCoord_; size_t api_yElemCoord_n_; vector2ptr(yElemCoord, &api_yElemCoord_, &api_yElemCoord_n_);
      double *api_zElemCoord_; size_t api_zElemCoord_n_; vector2ptr(zElemCoord, &api_zElemCoord_, &api_zElemCoord_n_);
      gmshViewProbe(tag, x, y, z, &api_value_, &api_value_n_, &distance, step, numComp, (int)gradient, distanceMax, api_xElemCoord_, api_xElemCoord_n_, api_yElemCoord_, api_yElemCoord_n_, api_zElemCoord_, api_zElemCoord_n_, dim, &ierr);
      if(ierr) throwLastError();
      value.assign(api_value_, api_value_ + api_value_n_); gmshFree(api_value_);
      gmshFree(api_xElemCoord_);
      gmshFree(api_yElemCoord_);
      gmshFree(api_zElemCoord_);
    }

    // Write the view to a file `fileName'. The export format is determined by the
    // file extension. Append to the file if `append' is set.
    inline void write(const int tag,
                      const std::string & fileName,
                      const bool append = false)
    {
      int ierr = 0;
      gmshViewWrite(tag, fileName.c_str(), (int)append, &ierr);
      if(ierr) throwLastError();
    }

    // Set the global visibility of the view `tag' per window to `value', where
    // `windowIndex' identifies the window in the window list.
    inline void setVisibilityPerWindow(const int tag,
                                       const int value,
                                       const int windowIndex = 0)
    {
      int ierr = 0;
      gmshViewSetVisibilityPerWindow(tag, value, windowIndex, &ierr);
      if(ierr) throwLastError();
    }

    namespace option { // View option handling functions

      // Set the numerical option `name' to value `value' for the view with tag
      // `tag'.
      inline void setNumber(const int tag,
                            const std::string & name,
                            const double value)
      {
        int ierr = 0;
        gmshViewOptionSetNumber(tag, name.c_str(), value, &ierr);
        if(ierr) throwLastError();
      }

      // Get the `value' of the numerical option `name' for the view with tag
      // `tag'.
      inline void getNumber(const int tag,
                            const std::string & name,
                            double & value)
      {
        int ierr = 0;
        gmshViewOptionGetNumber(tag, name.c_str(), &value, &ierr);
        if(ierr) throwLastError();
      }

      // Set the string option `name' to value `value' for the view with tag `tag'.
      inline void setString(const int tag,
                            const std::string & name,
                            const std::string & value)
      {
        int ierr = 0;
        gmshViewOptionSetString(tag, name.c_str(), value.c_str(), &ierr);
        if(ierr) throwLastError();
      }

      // Get the `value' of the string option `name' for the view with tag `tag'.
      inline void getString(const int tag,
                            const std::string & name,
                            std::string & value)
      {
        int ierr = 0;
        char *api_value_;
        gmshViewOptionGetString(tag, name.c_str(), &api_value_, &ierr);
        if(ierr) throwLastError();
        value = std::string(api_value_); gmshFree(api_value_);
      }

      // Set the color option `name' to the RGBA value (`r', `g', `b', `a') for the
      // view with tag `tag', where where `r', `g', `b' and `a' should be integers
      // between 0 and 255.
      inline void setColor(const int tag,
                           const std::string & name,
                           const int r,
                           const int g,
                           const int b,
                           const int a = 255)
      {
        int ierr = 0;
        gmshViewOptionSetColor(tag, name.c_str(), r, g, b, a, &ierr);
        if(ierr) throwLastError();
      }

      // Get the `r', `g', `b', `a' value of the color option `name' for the view
      // with tag `tag'.
      inline void getColor(const int tag,
                           const std::string & name,
                           int & r,
                           int & g,
                           int & b,
                           int & a)
      {
        int ierr = 0;
        gmshViewOptionGetColor(tag, name.c_str(), &r, &g, &b, &a, &ierr);
        if(ierr) throwLastError();
      }

      // Copy the options from the view with tag `refTag' to the view with tag
      // `tag'.
      inline void copy(const int refTag,
                       const int tag)
      {
        int ierr = 0;
        gmshViewOptionCopy(refTag, tag, &ierr);
        if(ierr) throwLastError();
      }

    } // namespace option

  } // namespace view

  namespace plugin { // Plugin functions

    // Set the numerical option `option' to the value `value' for plugin `name'.
    inline void setNumber(const std::string & name,
                          const std::string & option,
                          const double value)
    {
      int ierr = 0;
      gmshPluginSetNumber(name.c_str(), option.c_str(), value, &ierr);
      if(ierr) throwLastError();
    }

    // Set the string option `option' to the value `value' for plugin `name'.
    inline void setString(const std::string & name,
                          const std::string & option,
                          const std::string & value)
    {
      int ierr = 0;
      gmshPluginSetString(name.c_str(), option.c_str(), value.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Run the plugin `name'. Return the tag of the created view (if any).
    inline int run(const std::string & name)
    {
      int ierr = 0;
      int result_api_ = gmshPluginRun(name.c_str(), &ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

  } // namespace plugin

  namespace graphics { // Graphics functions

    // Draw all the OpenGL scenes.
    inline void draw()
    {
      int ierr = 0;
      gmshGraphicsDraw(&ierr);
      if(ierr) throwLastError();
    }

  } // namespace graphics

  namespace fltk { // FLTK graphical user interface functions

    // Create the FLTK graphical user interface. Can only be called in the main
    // thread.
    inline void initialize()
    {
      int ierr = 0;
      gmshFltkInitialize(&ierr);
      if(ierr) throwLastError();
    }

    // Wait at most `time' seconds for user interface events and return. If `time'
    // < 0, wait indefinitely. First automatically create the user interface if it
    // has not yet been initialized. Can only be called in the main thread.
    inline void wait(const double time = -1.)
    {
      int ierr = 0;
      gmshFltkWait(time, &ierr);
      if(ierr) throwLastError();
    }

    // Update the user interface (potentially creating new widgets and windows).
    // First automatically create the user interface if it has not yet been
    // initialized. Can only be called in the main thread: use `awake("update")' to
    // trigger an update of the user interface from another thread.
    inline void update()
    {
      int ierr = 0;
      gmshFltkUpdate(&ierr);
      if(ierr) throwLastError();
    }

    // Awake the main user interface thread and process pending events, and
    // optionally perform an action (currently the only `action' allowed is
    // "update").
    inline void awake(const std::string & action = "")
    {
      int ierr = 0;
      gmshFltkAwake(action.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Block the current thread until it can safely modify the user interface.
    inline void lock()
    {
      int ierr = 0;
      gmshFltkLock(&ierr);
      if(ierr) throwLastError();
    }

    // Release the lock that was set using lock.
    inline void unlock()
    {
      int ierr = 0;
      gmshFltkUnlock(&ierr);
      if(ierr) throwLastError();
    }

    // Run the event loop of the graphical user interface, i.e. repeatedly call
    // `wait()'. First automatically create the user interface if it has not yet
    // been initialized. Can only be called in the main thread.
    inline void run()
    {
      int ierr = 0;
      gmshFltkRun(&ierr);
      if(ierr) throwLastError();
    }

    // Check if the user interface is available (e.g. to detect if it has been
    // closed).
    inline int isAvailable()
    {
      int ierr = 0;
      int result_api_ = gmshFltkIsAvailable(&ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

    // Select entities in the user interface. If `dim' is >= 0, return only the
    // entities of the specified dimension (e.g. points if `dim' == 0).
    inline int selectEntities(gmsh::vectorpair & dimTags,
                              const int dim = -1)
    {
      int ierr = 0;
      int *api_dimTags_; size_t api_dimTags_n_;
      int result_api_ = gmshFltkSelectEntities(&api_dimTags_, &api_dimTags_n_, dim, &ierr);
      if(ierr) throwLastError();
      dimTags.resize(api_dimTags_n_ / 2); for(size_t i = 0; i < api_dimTags_n_ / 2; ++i){ dimTags[i].first = api_dimTags_[i * 2 + 0]; dimTags[i].second = api_dimTags_[i * 2 + 1]; } gmshFree(api_dimTags_);
      return result_api_;
    }

    // Select elements in the user interface.
    inline int selectElements(std::vector<std::size_t> & elementTags)
    {
      int ierr = 0;
      size_t *api_elementTags_; size_t api_elementTags_n_;
      int result_api_ = gmshFltkSelectElements(&api_elementTags_, &api_elementTags_n_, &ierr);
      if(ierr) throwLastError();
      elementTags.assign(api_elementTags_, api_elementTags_ + api_elementTags_n_); gmshFree(api_elementTags_);
      return result_api_;
    }

    // Select views in the user interface.
    inline int selectViews(std::vector<int> & viewTags)
    {
      int ierr = 0;
      int *api_viewTags_; size_t api_viewTags_n_;
      int result_api_ = gmshFltkSelectViews(&api_viewTags_, &api_viewTags_n_, &ierr);
      if(ierr) throwLastError();
      viewTags.assign(api_viewTags_, api_viewTags_ + api_viewTags_n_); gmshFree(api_viewTags_);
      return result_api_;
    }

    // Split the current window horizontally (if `how' = "h") or vertically (if
    // `how' = "v"), using ratio `ratio'. If `how' = "u", restore a single window.
    inline void splitCurrentWindow(const std::string & how = "v",
                                   const double ratio = 0.5)
    {
      int ierr = 0;
      gmshFltkSplitCurrentWindow(how.c_str(), ratio, &ierr);
      if(ierr) throwLastError();
    }

    // Set the current window by speficying its index (starting at 0) in the list
    // of all windows. When new windows are created by splits, new windows are
    // appended at the end of the list.
    inline void setCurrentWindow(const int windowIndex = 0)
    {
      int ierr = 0;
      gmshFltkSetCurrentWindow(windowIndex, &ierr);
      if(ierr) throwLastError();
    }

    // Set a status message in the current window. If `graphics' is set, display
    // the message inside the graphic window instead of the status bar.
    inline void setStatusMessage(const std::string & message,
                                 const bool graphics = false)
    {
      int ierr = 0;
      gmshFltkSetStatusMessage(message.c_str(), (int)graphics, &ierr);
      if(ierr) throwLastError();
    }

    // Show context window for the entity of dimension `dim' and tag `tag'.
    inline void showContextWindow(const int dim,
                                  const int tag)
    {
      int ierr = 0;
      gmshFltkShowContextWindow(dim, tag, &ierr);
      if(ierr) throwLastError();
    }

    // Open the `name' item in the menu tree.
    inline void openTreeItem(const std::string & name)
    {
      int ierr = 0;
      gmshFltkOpenTreeItem(name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Close the `name' item in the menu tree.
    inline void closeTreeItem(const std::string & name)
    {
      int ierr = 0;
      gmshFltkCloseTreeItem(name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

  } // namespace fltk

  namespace onelab { // ONELAB server functions

    // Set one or more parameters in the ONELAB database, encoded in `format'.
    inline void set(const std::string & data,
                    const std::string & format = "json")
    {
      int ierr = 0;
      gmshOnelabSet(data.c_str(), format.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Get all the parameters (or a single one if `name' is specified) from the
    // ONELAB database, encoded in `format'.
    inline void get(std::string & data,
                    const std::string & name = "",
                    const std::string & format = "json")
    {
      int ierr = 0;
      char *api_data_;
      gmshOnelabGet(&api_data_, name.c_str(), format.c_str(), &ierr);
      if(ierr) throwLastError();
      data = std::string(api_data_); gmshFree(api_data_);
    }

    // Get the names of the parameters in the ONELAB database matching the `search'
    // regular expression. If `search' is empty, return all the names.
    inline void getNames(std::vector<std::string> & names,
                         const std::string & search = "")
    {
      int ierr = 0;
      char **api_names_; size_t api_names_n_;
      gmshOnelabGetNames(&api_names_, &api_names_n_, search.c_str(), &ierr);
      if(ierr) throwLastError();
      names.resize(api_names_n_); for(size_t i = 0; i < api_names_n_; ++i){ names[i] = std::string(api_names_[i]); gmshFree(api_names_[i]); } gmshFree(api_names_);
    }

    // Set the value of the number parameter `name' in the ONELAB database. Create
    // the parameter if it does not exist; update the value if the parameter
    // exists.
    inline void setNumber(const std::string & name,
                          const std::vector<double> & value)
    {
      int ierr = 0;
      double *api_value_; size_t api_value_n_; vector2ptr(value, &api_value_, &api_value_n_);
      gmshOnelabSetNumber(name.c_str(), api_value_, api_value_n_, &ierr);
      if(ierr) throwLastError();
      gmshFree(api_value_);
    }

    // Set the value of the string parameter `name' in the ONELAB database. Create
    // the parameter if it does not exist; update the value if the parameter
    // exists.
    inline void setString(const std::string & name,
                          const std::vector<std::string> & value)
    {
      int ierr = 0;
      char **api_value_; size_t api_value_n_; vectorstring2charptrptr(value, &api_value_, &api_value_n_);
      gmshOnelabSetString(name.c_str(), api_value_, api_value_n_, &ierr);
      if(ierr) throwLastError();
      for(size_t i = 0; i < api_value_n_; ++i){ gmshFree(api_value_[i]); } gmshFree(api_value_);
    }

    // Get the value of the number parameter `name' from the ONELAB database.
    // Return an empty vector if the parameter does not exist.
    inline void getNumber(const std::string & name,
                          std::vector<double> & value)
    {
      int ierr = 0;
      double *api_value_; size_t api_value_n_;
      gmshOnelabGetNumber(name.c_str(), &api_value_, &api_value_n_, &ierr);
      if(ierr) throwLastError();
      value.assign(api_value_, api_value_ + api_value_n_); gmshFree(api_value_);
    }

    // Get the value of the string parameter `name' from the ONELAB database.
    // Return an empty vector if the parameter does not exist.
    inline void getString(const std::string & name,
                          std::vector<std::string> & value)
    {
      int ierr = 0;
      char **api_value_; size_t api_value_n_;
      gmshOnelabGetString(name.c_str(), &api_value_, &api_value_n_, &ierr);
      if(ierr) throwLastError();
      value.resize(api_value_n_); for(size_t i = 0; i < api_value_n_; ++i){ value[i] = std::string(api_value_[i]); gmshFree(api_value_[i]); } gmshFree(api_value_);
    }

    // Clear the ONELAB database, or remove a single parameter if `name' is given.
    inline void clear(const std::string & name = "")
    {
      int ierr = 0;
      gmshOnelabClear(name.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Run a ONELAB client. If `name' is provided, create a new ONELAB client with
    // name `name' and executes `command'. If not, try to run a client that might
    // be linked to the processed input files.
    inline void run(const std::string & name = "",
                    const std::string & command = "")
    {
      int ierr = 0;
      gmshOnelabRun(name.c_str(), command.c_str(), &ierr);
      if(ierr) throwLastError();
    }

  } // namespace onelab

  namespace logger { // Information logging functions

    // Write a `message'. `level' can be "info", "warning" or "error".
    inline void write(const std::string & message,
                      const std::string & level = "info")
    {
      int ierr = 0;
      gmshLoggerWrite(message.c_str(), level.c_str(), &ierr);
      if(ierr) throwLastError();
    }

    // Start logging messages.
    inline void start()
    {
      int ierr = 0;
      gmshLoggerStart(&ierr);
      if(ierr) throwLastError();
    }

    // Get logged messages.
    inline void get(std::vector<std::string> & log)
    {
      int ierr = 0;
      char **api_log_; size_t api_log_n_;
      gmshLoggerGet(&api_log_, &api_log_n_, &ierr);
      if(ierr) throwLastError();
      log.resize(api_log_n_); for(size_t i = 0; i < api_log_n_; ++i){ log[i] = std::string(api_log_[i]); gmshFree(api_log_[i]); } gmshFree(api_log_);
    }

    // Stop logging messages.
    inline void stop()
    {
      int ierr = 0;
      gmshLoggerStop(&ierr);
      if(ierr) throwLastError();
    }

    // Return wall clock time.
    inline double getWallTime()
    {
      int ierr = 0;
      double result_api_ = gmshLoggerGetWallTime(&ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

    // Return CPU time.
    inline double getCpuTime()
    {
      int ierr = 0;
      double result_api_ = gmshLoggerGetCpuTime(&ierr);
      if(ierr) throwLastError();
      return result_api_;
    }

    // Return last error message, if any.
    inline void getLastError(std::string & error)
    {
      int ierr = 0;
      char *api_error_;
      gmshLoggerGetLastError(&api_error_, &ierr);
      if(ierr) throw "Could not get last error";
      error = std::string(api_error_); gmshFree(api_error_);
    }

  } // namespace logger

} // namespace gmsh

#endif