Skip to content
Snippets Groups Projects
Select Git revision
  • f55633bad0503ecf844c3a10072aaf1d9e92e2f8
  • master default protected
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • alphashapes
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

GeoStringInterface.h

Blame
  • meshGRegionExtruded.cpp 20.64 KiB
    // Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
    //
    // See the LICENSE.txt file for license information. Please report all
    // issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
    
    #include <set>
    #include "GmshConfig.h"
    #include "GmshMessage.h"
    #include "GModel.h"
    #include "MTriangle.h"
    #include "MQuadrangle.h"
    #include "MTetrahedron.h"
    #include "MHexahedron.h"
    #include "MPrism.h"
    #include "MPyramid.h"
    #include "ExtrudeParams.h"
    #include "meshGFace.h"
    #include "meshGRegion.h"
    #include "Context.h"
    #include "MVertexRTree.h"
    
    #if defined(HAVE_QUADTRI)
    #include "QuadTriExtruded3D.h"
    #endif
    
    static void addTetrahedron(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
                               GRegion *to)
    {
      to->tetrahedra.push_back(new MTetrahedron(v1, v2, v3, v4));
    }
    
    static void addPyramid(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
                           MVertex *v5, GRegion *to)
    {
      to->pyramids.push_back(new MPyramid(v1, v2, v3, v4, v5));
    }
    
    static void addPrism(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
                         MVertex *v5, MVertex *v6, GRegion *to)
    {
      to->prisms.push_back(new MPrism(v1, v2, v3, v4, v5, v6));
    }
    
    static void addHexahedron(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
                              MVertex *v5, MVertex *v6, MVertex *v7, MVertex *v8,
                              GRegion *to)
    {
      to->hexahedra.push_back(new MHexahedron(v1, v2, v3, v4, v5, v6, v7, v8));
    }
    
    static void createPriPyrTet(std::vector<MVertex *> &v, GRegion *to,
                                MElement *source)
    {
      int dup[3];
      int j = 0;
      for(int i = 0; i < 3; i++)
        if(v[i] == v[i + 3]) dup[j++] = i;
    
      if(j == 2) {
        if(dup[0] == 0 && dup[1] == 1)
          addTetrahedron(v[0], v[1], v[2], v[5], to);
        else if(dup[0] == 1 && dup[1] == 2)
          addTetrahedron(v[0], v[1], v[2], v[3], to);
        else
          addTetrahedron(v[0], v[1], v[2], v[4], to);
      }
      else if(j == 1) {
        if(dup[0] == 0)
          addPyramid(v[1], v[4], v[5], v[2], v[0], to);
        else if(dup[0] == 1)
          addPyramid(v[0], v[2], v[5], v[3], v[1], to);
        else
          addPyramid(v[0], v[1], v[4], v[3], v[2], to);
      }
      else {
        addPrism(v[0], v[1], v[2], v[3], v[4], v[5], to);
        if(j) Msg::Error("Degenerated prism in extrusion of volume %d", to->tag());
      }
    }
    
    static void createHexPri(std::vector<MVertex *> &v, GRegion *to,
                             MElement *source)
    {
      int dup[4];
      int j = 0;
      for(int i = 0; i < 4; i++)
        if(v[i] == v[i + 4]) dup[j++] = i;
    
      if(j == 2) {
        if(dup[0] == 0 && dup[1] == 1)
          addPrism(v[0], v[3], v[7], v[1], v[2], v[6], to);
        else if(dup[0] == 1 && dup[1] == 2)
          addPrism(v[0], v[1], v[4], v[3], v[2], v[7], to);
        else if(dup[0] == 2 && dup[1] == 3)
          addPrism(v[0], v[3], v[4], v[1], v[2], v[5], to);
        else if(dup[0] == 0 && dup[1] == 3)
          addPrism(v[0], v[1], v[5], v[3], v[2], v[6], to);
        else
          Msg::Error("Uncoherent hexahedron in extrusion of volume %d", to->tag());
      }
      else {
        addHexahedron(v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7], to);
        if(j)
          Msg::Error("Degenerated hexahedron in extrusion of volume %d", to->tag());
      }
    }
    
    static void createTet(MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
                          GRegion *to, MElement *source)
    {
      if(v1 != v2 && v1 != v3 && v1 != v4 && v2 != v3 && v2 != v4 && v3 != v4)
        addTetrahedron(v1, v2, v3, v4, to);
    }
    
    static int getExtrudedVertices(MElement *ele, ExtrudeParams *ep, int j, int k,
                                   MVertexRTree &pos, std::vector<MVertex *> &verts)
    {
      double x[8], y[8], z[8];
      int n = ele->getNumVertices();
      for(int p = 0; p < n; p++) {
        MVertex *v = ele->getVertex(p);
        x[p] = x[p + n] = v->x();
        y[p] = y[p + n] = v->y();
        z[p] = z[p + n] = v->z();
      }
      for(int p = 0; p < n; p++) {
        ep->Extrude(j, k, x[p], y[p], z[p]);
        ep->Extrude(j, k + 1, x[p + n], y[p + n], z[p + n]);
      }
      for(int p = 0; p < 2 * n; p++) {
        MVertex *tmp = pos.find(x[p], y[p], z[p]);
        if(!tmp)
          Msg::Error("Could not find extruded vertex (%.16g, %.16g, %.16g)", x[p],
                     y[p], z[p]);
        else
          verts.push_back(tmp);
      }
      return verts.size();
    }
    
    static void extrudeMesh(GFace *from, GRegion *to, MVertexRTree &pos)
    {
      ExtrudeParams *ep = to->meshAttributes.extrude;
    
      // interior vertices
      std::vector<MVertex *> mesh_vertices = from->mesh_vertices;
    
      // add all embedded vertices
      std::vector<MVertex *> embedded = from->getEmbeddedMeshVertices();
      mesh_vertices.insert(mesh_vertices.end(), embedded.begin(), embedded.end());
    
      // add all vertices on surface seams
      std::set<MVertex *> seam;
      std::vector<GEdge *> const &l_edges = from->edges();
      for(std::vector<GEdge *>::const_iterator it = l_edges.begin();
          it != l_edges.end(); ++it) {
        if((*it)->isSeam(from)) {
          seam.insert((*it)->mesh_vertices.begin(), (*it)->mesh_vertices.end());
          if((*it)->getBeginVertex())
            seam.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
                        (*it)->getBeginVertex()->mesh_vertices.end());
          if((*it)->getEndVertex())
            seam.insert((*it)->getEndVertex()->mesh_vertices.begin(),
                        (*it)->getEndVertex()->mesh_vertices.end());
        }
      }
      mesh_vertices.insert(mesh_vertices.end(), seam.begin(), seam.end());
    
      // create extruded vertices
      for(std::size_t i = 0; i < mesh_vertices.size(); i++) {
        MVertex *v = mesh_vertices[i];
        for(int j = 0; j < ep->mesh.NbLayer; j++) {
          for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
            double x = v->x(), y = v->y(), z = v->z();
            ep->Extrude(j, k + 1, x, y, z);
            if(j != ep->mesh.NbLayer - 1 || k != ep->mesh.NbElmLayer[j] - 1) {
              MVertex *newv = new MVertex(x, y, z, to);
              to->mesh_vertices.push_back(newv);
              pos.insert(newv);
            }
          }
        }
      }
    
    #if defined(HAVE_QUADTRI)
      if(ep && ep->mesh.ExtrudeMesh && ep->mesh.QuadToTri && ep->mesh.Recombine) {
        meshQuadToTriRegion(to, pos);
        return;
      }
    #endif
    
      // create elements (note that it would be faster to access the *interior*
      // nodes by direct indexing, but it's just simpler to query everything by
      // position)
      for(std::size_t i = 0; i < from->triangles.size(); i++) {
        for(int j = 0; j < ep->mesh.NbLayer; j++) {
          for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
            std::vector<MVertex *> verts;
            if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, verts) == 6) {
              createPriPyrTet(verts, to, from->triangles[i]);
            }
          }
        }
      }
    
      if(from->quadrangles.size() && !ep->mesh.Recombine) {
        Msg::Error("Cannot extrude quadrangles without Recombine");
      }
      else {
        for(std::size_t i = 0; i < from->quadrangles.size(); i++) {
          for(int j = 0; j < ep->mesh.NbLayer; j++) {
            for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
              std::vector<MVertex *> verts;
              if(getExtrudedVertices(from->quadrangles[i], ep, j, k, pos, verts) ==
                 8)
                createHexPri(verts, to, from->quadrangles[i]);
            }
          }
        }
      }
    }
    
    static void insertAllVertices(GRegion *gr, MVertexRTree &pos)
    {
      pos.insert(gr->mesh_vertices);
      std::vector<GFace *> faces = gr->faces();
      std::vector<GFace *>::iterator itf = faces.begin();
      while(itf != faces.end()) {
        pos.insert((*itf)->mesh_vertices);
        std::vector<MVertex *> embedded = (*itf)->getEmbeddedMeshVertices();
        pos.insert(embedded);
        std::vector<GEdge *> const &edges = (*itf)->edges();
        std::vector<GEdge *>::const_iterator ite = edges.begin();
        while(ite != edges.end()) {
          pos.insert((*ite)->mesh_vertices);
          if((*ite)->getBeginVertex())
            pos.insert((*ite)->getBeginVertex()->mesh_vertices);
          if((*ite)->getEndVertex())
             pos.insert((*ite)->getEndVertex()->mesh_vertices);
          ++ite;
        }
        ++itf;
      }
    }
    
    void meshGRegionExtruded::operator()(GRegion *gr)
    {
      gr->model()->setCurrentMeshEntity(gr);
    
      if(gr->geomType() == GEntity::DiscreteVolume) return;
    
      ExtrudeParams *ep = gr->meshAttributes.extrude;
    
      if(!ep || !ep->mesh.ExtrudeMesh || ep->geo.Mode != EXTRUDED_ENTITY) return;
    
      Msg::Info("Meshing volume %d (extruded)", gr->tag());
    
      // destroy the mesh if it exists
      deMeshGRegion dem;
      dem(gr);
    
      // build an rtree with all the vertices on the boundary of gr
      MVertexRTree pos(CTX::instance()->geom.tolerance * CTX::instance()->lc);
      insertAllVertices(gr, pos);
    
      // volume is extruded from a surface
      GFace *from = gr->model()->getFaceByTag(std::abs(ep->geo.Source));
      if(!from) {
        Msg::Error("Unknown source surface %d for extrusion", ep->geo.Source);
        return;
      }
    
      extrudeMesh(from, gr, pos);
    
      // carve holes if any (only do it now if the mesh is final, i.e., if
      // the mesh is recombined)
      if(ep->mesh.Holes.size() && ep->mesh.Recombine) {
        std::map<int, std::pair<double, std::vector<int> > >::iterator it;
        for(it = ep->mesh.Holes.begin(); it != ep->mesh.Holes.end(); it++)
          carveHole(gr, it->first, it->second.first, it->second.second);
      }
    }
    
    static int edgeExists(MVertex *v1, MVertex *v2,
                          std::set<std::pair<MVertex *, MVertex *> > &edges)
    {
      std::pair<MVertex *, MVertex *> p(std::min(v1, v2), std::max(v1, v2));
      return edges.count(p);
    }
    
    static void createEdge(MVertex *v1, MVertex *v2,
                           std::set<std::pair<MVertex *, MVertex *> > &edges)
    {
      std::pair<MVertex *, MVertex *> p(std::min(v1, v2), std::max(v1, v2));
      edges.insert(p);
    }
    
    static void deleteEdge(MVertex *v1, MVertex *v2,
                           std::set<std::pair<MVertex *, MVertex *> > &edges)
    {
      std::pair<MVertex *, MVertex *> p(std::min(v1, v2), std::max(v1, v2));
      edges.erase(p);
    }
    
    // subdivide the 3 lateral faces of each prism
    static void phase1(GRegion *gr, MVertexRTree &pos,
                       std::set<std::pair<MVertex *, MVertex *> > &edges)
    {
      ExtrudeParams *ep = gr->meshAttributes.extrude;
      GFace *from = gr->model()->getFaceByTag(std::abs(ep->geo.Source));
      if(!from) return;
    
      for(std::size_t i = 0; i < from->triangles.size(); i++) {
        for(int j = 0; j < ep->mesh.NbLayer; j++) {
          for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
            std::vector<MVertex *> v;
            if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, v) == 6) {
    #if 0 // old
              if(!edgeExists(v[0], v[4], edges))
                createEdge(v[1], v[3], edges);
              if(!edgeExists(v[4], v[2], edges))
                createEdge(v[1], v[5], edges);
              if(!edgeExists(v[3], v[2], edges))
                createEdge(v[0], v[5], edges);
    #else // new from Michel Benhamou
              if(v[1] < v[0])
                createEdge(v[1], v[3], edges);
              else
                createEdge(v[0], v[4], edges);
              if(v[1] < v[2])
                createEdge(v[1], v[5], edges);
              else
                createEdge(v[4], v[2], edges);
              if(v[0] < v[2])
                createEdge(v[0], v[5], edges);
              else
                createEdge(v[3], v[2], edges);
    #endif
            }
          }
        }
      }
    }
    
    // modify lateral edges to make them "tet-compatible"
    static void phase2(GRegion *gr, MVertexRTree &pos,
                       std::set<std::pair<MVertex *, MVertex *> > &edges,
                       std::set<std::pair<MVertex *, MVertex *> > &edges_swap,
                       int &swap)
    {
      ExtrudeParams *ep = gr->meshAttributes.extrude;
      GFace *from = gr->model()->getFaceByTag(std::abs(ep->geo.Source));
      if(!from) return;
    
      for(std::size_t i = 0; i < from->triangles.size(); i++) {
        for(int j = 0; j < ep->mesh.NbLayer; j++) {
          for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
            std::vector<MVertex *> v;
            if(getExtrudedVertices(from->triangles[i], ep, j, k, pos, v) == 6) {
              if(edgeExists(v[3], v[1], edges) && edgeExists(v[4], v[2], edges) &&
                 edgeExists(v[0], v[5], edges)) {
                swap++;
                if(!edgeExists(v[3], v[1], edges_swap)) {
                  deleteEdge(v[3], v[1], edges);
                  createEdge(v[0], v[4], edges);
                  createEdge(v[3], v[1], edges_swap);
                  createEdge(v[0], v[4], edges_swap);
                }
                else if(!edgeExists(v[4], v[2], edges_swap)) {
                  deleteEdge(v[4], v[2], edges);
                  createEdge(v[1], v[5], edges);
                  createEdge(v[4], v[2], edges_swap);
                  createEdge(v[1], v[5], edges_swap);
                }
                else if(!edgeExists(v[0], v[5], edges_swap)) {
                  deleteEdge(v[0], v[5], edges);
                  createEdge(v[3], v[2], edges);
                  createEdge(v[0], v[5], edges_swap);
                  createEdge(v[3], v[2], edges_swap);
                }
              }
              else if(edgeExists(v[0], v[4], edges) &&
                      edgeExists(v[1], v[5], edges) &&
                      edgeExists(v[3], v[2], edges)) {
                swap++;
                if(!edgeExists(v[0], v[4], edges_swap)) {
                  deleteEdge(v[0], v[4], edges);
                  createEdge(v[3], v[1], edges);
                  createEdge(v[0], v[4], edges_swap);
                  createEdge(v[3], v[1], edges_swap);
                }
                else if(!edgeExists(v[1], v[5], edges_swap)) {
                  deleteEdge(v[1], v[5], edges);
                  createEdge(v[4], v[2], edges);
                  createEdge(v[1], v[5], edges_swap);
                  createEdge(v[4], v[2], edges_swap);
                }
                else if(!edgeExists(v[3], v[2], edges_swap)) {
                  deleteEdge(v[3], v[2], edges);
                  createEdge(v[0], v[5], edges);
                  createEdge(v[3], v[2], edges_swap);
                  createEdge(v[0], v[5], edges_swap);
                }
              }
            }
          }
        }
      }
    }
    
    // create tets
    static void phase3(GRegion *gr, MVertexRTree &pos,
                       std::set<std::pair<MVertex *, MVertex *> > &edges)
    {
      ExtrudeParams *ep = gr->meshAttributes.extrude;
      GFace *from = gr->model()->getFaceByTag(std::abs(ep->geo.Source));
      if(!from) return;
    
      for(std::size_t i = 0; i < from->triangles.size(); i++) {
        MTriangle *tri = from->triangles[i];
        for(int j = 0; j < ep->mesh.NbLayer; j++) {
          for(int k = 0; k < ep->mesh.NbElmLayer[j]; k++) {
            std::vector<MVertex *> v;
            if(getExtrudedVertices(tri, ep, j, k, pos, v) == 6) {
              if(edgeExists(v[3], v[1], edges) && edgeExists(v[4], v[2], edges) &&
                 edgeExists(v[3], v[2], edges)) {
                createTet(v[0], v[1], v[2], v[3], gr, tri);
                createTet(v[3], v[4], v[5], v[2], gr, tri);
                createTet(v[1], v[3], v[4], v[2], gr, tri);
              }
              else if(edgeExists(v[3], v[1], edges) &&
                      edgeExists(v[1], v[5], edges) &&
                      edgeExists(v[3], v[2], edges)) {
                createTet(v[0], v[1], v[2], v[3], gr, tri);
                createTet(v[3], v[4], v[5], v[1], gr, tri);
                createTet(v[3], v[1], v[5], v[2], gr, tri);
              }
              else if(edgeExists(v[3], v[1], edges) &&
                      edgeExists(v[1], v[5], edges) &&
                      edgeExists(v[5], v[0], edges)) {
                createTet(v[0], v[1], v[2], v[5], gr, tri);
                createTet(v[3], v[4], v[5], v[1], gr, tri);
                createTet(v[1], v[3], v[5], v[0], gr, tri);
              }
              else if(edgeExists(v[4], v[0], edges) &&
                      edgeExists(v[4], v[2], edges) &&
                      edgeExists(v[3], v[2], edges)) {
                createTet(v[0], v[1], v[2], v[4], gr, tri);
                createTet(v[3], v[4], v[5], v[2], gr, tri);
                createTet(v[0], v[3], v[4], v[2], gr, tri);
              }
              else if(edgeExists(v[4], v[0], edges) &&
                      edgeExists(v[4], v[2], edges) &&
                      edgeExists(v[5], v[0], edges)) {
                createTet(v[0], v[1], v[2], v[4], gr, tri);
                createTet(v[3], v[4], v[5], v[0], gr, tri);
                createTet(v[0], v[2], v[4], v[5], gr, tri);
              }
              else if(edgeExists(v[4], v[0], edges) &&
                      edgeExists(v[1], v[5], edges) &&
                      edgeExists(v[5], v[0], edges)) {
                createTet(v[0], v[1], v[2], v[5], gr, tri);
                createTet(v[3], v[4], v[5], v[0], gr, tri);
                createTet(v[0], v[1], v[4], v[5], gr, tri);
              }
            }
          }
        }
      }
    }
    
    int SubdivideExtrudedMesh(GModel *m)
    {
      // get all non-recombined extruded regions and vertices; also, create a vector
      // of quadToTri regions that have NOT been meshed yet
      std::vector<GRegion *> regions;
    #if defined(HAVE_QUADTRI)
      std::vector<GRegion *> regions_quadToTri;
    #endif
    
      MVertexRTree pos(CTX::instance()->geom.tolerance * CTX::instance()->lc);
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); it++) {
        ExtrudeParams *ep = (*it)->meshAttributes.extrude;
        if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY &&
           !ep->mesh.Recombine) {
          regions.push_back(*it);
          insertAllVertices(*it, pos);
        }
    #if defined(HAVE_QUADTRI)
        // create vector of valid quadToTri regions...not all will necessarily be
        // meshed here.
        if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY &&
           ep->mesh.Recombine && ep->mesh.QuadToTri) {
          regions_quadToTri.push_back(*it);
        }
    #endif
      }
    
      if(regions.empty()) return 0;
      Msg::Info("Subdividing extruded mesh");
    
      // create edges on lateral sides of "prisms"
      std::set<std::pair<MVertex *, MVertex *> > edges;
      for(std::size_t i = 0; i < regions.size(); i++)
        phase1(regions[i], pos, edges);
    
      // swap lateral edges to make them "tet-compatible"
      int j = 0, swap;
      std::set<std::pair<MVertex *, MVertex *> > edges_swap;
      do {
        swap = 0;
        for(std::size_t i = 0; i < regions.size(); i++)
          phase2(regions[i], pos, edges, edges_swap, swap);
        Msg::Info("Swapping %d", swap);
        if(j && j == swap) {
          Msg::Error("Unable to subdivide extruded mesh: change surface mesh or");
          Msg::Error("recombine extrusion instead");
          return -1;
        }
        j = swap;
      } while(swap);
    
      // delete volume elements and create tetrahedra instead
      for(std::size_t i = 0; i < regions.size(); i++) {
        GRegion *gr = regions[i];
    
        for(std::size_t i = 0; i < gr->tetrahedra.size(); i++)
          delete gr->tetrahedra[i];
        gr->tetrahedra.clear();
        for(std::size_t i = 0; i < gr->hexahedra.size(); i++)
          delete gr->hexahedra[i];
        gr->hexahedra.clear();
        for(std::size_t i = 0; i < gr->prisms.size(); i++) delete gr->prisms[i];
        gr->prisms.clear();
        for(std::size_t i = 0; i < gr->pyramids.size(); i++)
          delete gr->pyramids[i];
        gr->pyramids.clear();
        phase3(gr, pos, edges);
    
        // re-Extrude bounding surfaces using edges as constraint
        std::vector<GFace *> faces = gr->faces();
        for(std::vector<GFace *>::iterator it = faces.begin(); it != faces.end();
            it++) {
          ExtrudeParams *ep = (*it)->meshAttributes.extrude;
          if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == EXTRUDED_ENTITY &&
             !ep->mesh.Recombine) {
            GFace *gf = *it;
            Msg::Info("Remeshing surface %d", gf->tag());
            for(std::size_t i = 0; i < gf->triangles.size(); i++)
              delete gf->triangles[i];
            gf->triangles.clear();
            for(std::size_t i = 0; i < gf->quadrangles.size(); i++)
              delete gf->quadrangles[i];
            gf->quadrangles.clear();
            MeshExtrudedSurface(gf, &edges);
          }
        }
      }
    
    #if defined(HAVE_QUADTRI)
      // now mesh the QuadToTri regions. Everything can be done locally for each
      // quadToTri region, but still use edge set from above just to make sure
      // laterals get remeshed properly ( QuadToTriEdgeGenerator detects if the
      // neighbor has been meshed or if a lateral surface should remain static for
      // any other reason).  If this function detects allNonGlobalSharedLaterals, it
      // won't mesh the region (should already be done in ExtrudeMesh).
      for(std::size_t i = 0; i < regions_quadToTri.size(); i++) {
        GRegion *gr = regions_quadToTri[i];
        MVertexRTree pos_local(CTX::instance()->geom.tolerance *
                               CTX::instance()->lc);
        insertAllVertices(gr, pos_local);
        meshQuadToTriRegionAfterGlobalSubdivide(gr, &edges, pos_local);
      }
    #endif
    
      // carve holes if any (TODO: update extrusion information)
      for(std::size_t i = 0; i < regions.size(); i++) {
        GRegion *gr = regions[i];
        ExtrudeParams *ep = gr->meshAttributes.extrude;
        if(ep->mesh.Holes.size()) {
          std::map<int, std::pair<double, std::vector<int> > >::iterator it;
          for(it = ep->mesh.Holes.begin(); it != ep->mesh.Holes.end(); it++)
            carveHole(gr, it->first, it->second.first, it->second.second);
        }
      }
    
    #if defined(HAVE_QUADTRI)
      for(std::size_t i = 0; i < regions_quadToTri.size(); i++) {
        GRegion *gr = regions_quadToTri[i];
        ExtrudeParams *ep = gr->meshAttributes.extrude;
        if(ep->mesh.Holes.size()) {
          std::map<int, std::pair<double, std::vector<int> > >::iterator it;
          for(it = ep->mesh.Holes.begin(); it != ep->mesh.Holes.end(); it++)
            carveHole(gr, it->first, it->second.first, it->second.second);
        }
      }
    #endif
    
      return 1;
    }