Skip to content
Snippets Groups Projects
Select Git revision
  • 83a1a852b7405612e563efb02171b0ceb3bc2f5e
  • master default protected
  • overlaps_tags_and_distributed_export
  • overlaps_tags_and_distributed_export_rebased
  • relaying
  • alphashapes
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 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

InnerVertexPlacement.cpp

Blame
  • InnerVertexPlacement.cpp 22.32 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.
    //
    // Contributed by Amaury Johnen
    
    #include "InnerVertexPlacement.h"
    #include "GmshDefines.h"
    #include "pointsGenerators.h"
    
    namespace {
      typedef std::pair<std::pair<int, int>, int> mytuple;
    
      mytuple make_mytuple(int a, int b, int c)
      {
        return std::make_pair(std::make_pair(a, b), c);
      }
    
      void double2int(const fullMatrix<double> &doubleMat, fullMatrix<int> &intMat)
      {
        intMat.resize(doubleMat.size1(), doubleMat.size2(), false);
        for(int i = 0; i < intMat.size1(); ++i) {
          for(int j = 0; j < intMat.size2(); ++j) {
            intMat(i, j) = static_cast<int>(doubleMat(i, j) + .5);
          }
        }
      }
    
      std::map<int, fullMatrix<double> *> storedMatrices[6];
    } // namespace
    
    fullMatrix<double> *getInnerVertexPlacement(int type, int order)
    {
      if(type < 3 || type > 8) return NULL;
      std::map<int, fullMatrix<double> *>::iterator it;
      it = storedMatrices[type - 3].find(order);
      if(it != storedMatrices[type - 3].end()) {
        return it->second;
      }
      else {
        fullMatrix<double> *matrix = new fullMatrix<double>();
        switch(type) {
        case TYPE_TRI:
          *matrix = gmshGenerateInnerVertexPlacementTriangle(order);
          break;
        case TYPE_QUA:
          *matrix = gmshGenerateInnerVertexPlacementQuadrangle(order);
          break;
        case TYPE_TET:
          *matrix = gmshGenerateInnerVertexPlacementTetrahedron(order);
          break;
        case TYPE_PRI:
          *matrix = gmshGenerateInnerVertexPlacementPrism(order);
          break;
        case TYPE_HEX:
          *matrix = gmshGenerateInnerVertexPlacementHexahedron(order);
          break;
        case TYPE_PYR:
          *matrix = gmshGenerateInnerVertexPlacementPyramid(order);
          break;
        }
        storedMatrices[type - 3][order] = matrix;
        return matrix;
      }
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementTriangle(int order)
    {
      if(order < 3) return fullMatrix<double>(0, 0);
    
      const int szInc = 3 * order;
      const int szComp = (order + 1) * (order + 2) / 2;
      const int szInt = szComp - szInc;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsTriangle(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<std::pair<int, int>, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        coord2idx[std::make_pair(u, v)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        int q = n - u - v;
        double xi = (double)u / n;
        double eta = (double)v / n;
        double rho = (double)q / n;
    
        M(i, coord2idx[std::make_pair(u + v, 0)]) += xi;
        M(i, coord2idx[std::make_pair(u + q, v)]) += xi;
        M(i, coord2idx[std::make_pair(0, v + u)]) += eta;
        M(i, coord2idx[std::make_pair(u, v + q)]) += eta;
        M(i, coord2idx[std::make_pair(0, v)]) += rho;
        M(i, coord2idx[std::make_pair(u, 0)]) += rho;
        M(i, coord2idx[std::make_pair(n, 0)]) -= xi;
        M(i, coord2idx[std::make_pair(0, n)]) -= eta;
        M(i, coord2idx[std::make_pair(0, 0)]) -= rho;
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementQuadrangle(int order)
    {
      if(order < 2) return fullMatrix<double>(0, 0);
    
      const int szInc = 4 * order;
      const int szComp = (order + 1) * (order + 1);
      const int szInt = szComp - szInc;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsQuadrangle(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<std::pair<int, int>, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        coord2idx[std::make_pair(u, v)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        double xi = (double)u / n;
        double eta = (double)v / n;
    
        M(i, coord2idx[std::make_pair(0, v)]) += 1 - xi;
        M(i, coord2idx[std::make_pair(n, v)]) += xi;
        M(i, coord2idx[std::make_pair(u, 0)]) += 1 - eta;
        M(i, coord2idx[std::make_pair(u, n)]) += eta;
        M(i, coord2idx[std::make_pair(0, 0)]) -= (1 - xi) * (1 - eta);
        M(i, coord2idx[std::make_pair(n, 0)]) -= xi * (1 - eta);
        M(i, coord2idx[std::make_pair(n, n)]) -= xi * eta;
        M(i, coord2idx[std::make_pair(0, n)]) -= (1 - xi) * eta;
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementTetrahedron(int order)
    {
      if(order < 4) return fullMatrix<double>(0, 0);
    
      const int szInt = (order - 3) * (order - 2) * (order - 1) / 6;
      const int szComp = (order + 1) * (order + 2) * (order + 3) / 6;
      const int szInc = szComp - szInt;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsTetrahedron(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<mytuple, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        int w = coordinates(i, 2);
        coord2idx[make_mytuple(u, v, w)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        int w = coordinates(szInc + i, 2);
        int q = n - u - v - w;
        double xi = (double)u / n;
        double eta = (double)v / n;
        double zeta = (double)w / n;
        double rho = (double)q / n;
    
        M(i, coord2idx[make_mytuple(u + v, 0, w)]) += xi;
        M(i, coord2idx[make_mytuple(u + w, v, 0)]) += xi;
        M(i, coord2idx[make_mytuple(u + q, v, w)]) += xi;
        M(i, coord2idx[make_mytuple(0, v + u, w)]) += eta;
        M(i, coord2idx[make_mytuple(u, v + w, 0)]) += eta;
        M(i, coord2idx[make_mytuple(u, v + q, w)]) += eta;
        M(i, coord2idx[make_mytuple(0, v, w + u)]) += zeta;
        M(i, coord2idx[make_mytuple(u, 0, w + v)]) += zeta;
        M(i, coord2idx[make_mytuple(u, v, w + q)]) += zeta;
        M(i, coord2idx[make_mytuple(0, v, w)]) += rho;
        M(i, coord2idx[make_mytuple(u, 0, w)]) += rho;
        M(i, coord2idx[make_mytuple(u, v, 0)]) += rho;
    
        M(i, coord2idx[make_mytuple(n - q, 0, 0)]) -= xi;
        M(i, coord2idx[make_mytuple(n - w, 0, w)]) -= xi;
        M(i, coord2idx[make_mytuple(n - v, v, 0)]) -= xi;
        M(i, coord2idx[make_mytuple(0, n - q, 0)]) -= eta;
        M(i, coord2idx[make_mytuple(0, n - w, w)]) -= eta;
        M(i, coord2idx[make_mytuple(u, n - u, 0)]) -= eta;
        M(i, coord2idx[make_mytuple(0, 0, n - q)]) -= zeta;
        M(i, coord2idx[make_mytuple(0, v, n - v)]) -= zeta;
        M(i, coord2idx[make_mytuple(u, 0, n - u)]) -= zeta;
        M(i, coord2idx[make_mytuple(0, 0, w)]) -= rho;
        M(i, coord2idx[make_mytuple(0, v, 0)]) -= rho;
        M(i, coord2idx[make_mytuple(u, 0, 0)]) -= rho;
    
        M(i, coord2idx[make_mytuple(n, 0, 0)]) += xi;
        M(i, coord2idx[make_mytuple(0, n, 0)]) += eta;
        M(i, coord2idx[make_mytuple(0, 0, n)]) += zeta;
        M(i, coord2idx[make_mytuple(0, 0, 0)]) += rho;
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementHexahedron(int order)
    {
      if(order < 2) return fullMatrix<double>(0, 0);
    
      const int szInt = (order - 1) * (order - 1) * (order - 1);
      const int szComp = (order + 1) * (order + 1) * (order + 1);
      const int szInc = szComp - szInt;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsHexahedron(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<mytuple, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        int w = coordinates(i, 2);
        coord2idx[make_mytuple(u, v, w)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        int w = coordinates(szInc + i, 2);
        double xi = (double)u / n;
        double eta = (double)v / n;
        double zeta = (double)w / n;
    
        M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - xi;
        M(i, coord2idx[make_mytuple(n, v, w)]) += xi;
        M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - eta;
        M(i, coord2idx[make_mytuple(u, n, w)]) += eta;
        M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - zeta;
        M(i, coord2idx[make_mytuple(u, v, n)]) += zeta;
    
        M(i, coord2idx[make_mytuple(0, 0, w)]) -= (1 - xi) * (1 - eta);
        M(i, coord2idx[make_mytuple(0, n, w)]) -= (1 - xi) * eta;
        M(i, coord2idx[make_mytuple(n, 0, w)]) -= xi * (1 - eta);
        M(i, coord2idx[make_mytuple(n, n, w)]) -= xi * eta;
        M(i, coord2idx[make_mytuple(u, 0, 0)]) -= (1 - eta) * (1 - zeta);
        M(i, coord2idx[make_mytuple(u, 0, n)]) -= (1 - eta) * zeta;
        M(i, coord2idx[make_mytuple(u, n, 0)]) -= eta * (1 - zeta);
        M(i, coord2idx[make_mytuple(u, n, n)]) -= eta * zeta;
        M(i, coord2idx[make_mytuple(0, v, 0)]) -= (1 - zeta) * (1 - xi);
        M(i, coord2idx[make_mytuple(n, v, 0)]) -= (1 - zeta) * xi;
        M(i, coord2idx[make_mytuple(0, v, n)]) -= zeta * (1 - xi);
        M(i, coord2idx[make_mytuple(n, v, n)]) -= zeta * xi;
    
        M(i, coord2idx[make_mytuple(0, 0, 0)]) += (1 - xi) * (1 - eta) * (1 - zeta);
        M(i, coord2idx[make_mytuple(0, 0, n)]) += (1 - xi) * (1 - eta) * zeta;
        M(i, coord2idx[make_mytuple(0, n, 0)]) += (1 - xi) * eta * (1 - zeta);
        M(i, coord2idx[make_mytuple(0, n, n)]) += (1 - xi) * eta * zeta;
        M(i, coord2idx[make_mytuple(n, 0, 0)]) += xi * (1 - eta) * (1 - zeta);
        M(i, coord2idx[make_mytuple(n, 0, n)]) += xi * (1 - eta) * zeta;
        M(i, coord2idx[make_mytuple(n, n, 0)]) += xi * eta * (1 - zeta);
        M(i, coord2idx[make_mytuple(n, n, n)]) += xi * eta * zeta;
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementPrism(int order)
    {
      if(order < 3) return fullMatrix<double>(0, 0);
    
      const int szInt = (order - 1) * (order - 2) * (order - 1) / 2;
      const int szComp = (order + 1) * (order + 1) * (order + 2) / 2;
      const int szInc = szComp - szInt;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsPrism(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<mytuple, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        int w = coordinates(i, 2);
        coord2idx[make_mytuple(u, v, w)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        int w = coordinates(szInc + i, 2);
        int q = n - u - v;
        double xi = (double)u / n;
        double eta = (double)v / n;
        double zeta = (double)w / n;
        double rho = (double)q / n;
    
        M(i, coord2idx[make_mytuple(u + v, 0, w)]) += xi;
        M(i, coord2idx[make_mytuple(u + q, v, w)]) += xi;
        M(i, coord2idx[make_mytuple(0, v + u, w)]) += eta;
        M(i, coord2idx[make_mytuple(u, v + q, w)]) += eta;
        M(i, coord2idx[make_mytuple(0, v, w)]) += rho;
        M(i, coord2idx[make_mytuple(u, 0, w)]) += rho;
        M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - zeta;
        M(i, coord2idx[make_mytuple(u, v, n)]) += zeta;
    
        M(i, coord2idx[make_mytuple(n, 0, w)]) -= xi;
        M(i, coord2idx[make_mytuple(0, n, w)]) -= eta;
        M(i, coord2idx[make_mytuple(0, 0, w)]) -= rho;
        M(i, coord2idx[make_mytuple(u + v, 0, 0)]) -= xi * (1 - zeta);
        M(i, coord2idx[make_mytuple(u + q, v, 0)]) -= xi * (1 - zeta);
        M(i, coord2idx[make_mytuple(0, v + u, 0)]) -= eta * (1 - zeta);
        M(i, coord2idx[make_mytuple(u, v + q, 0)]) -= eta * (1 - zeta);
        M(i, coord2idx[make_mytuple(0, v, 0)]) -= rho * (1 - zeta);
        M(i, coord2idx[make_mytuple(u, 0, 0)]) -= rho * (1 - zeta);
        M(i, coord2idx[make_mytuple(u + v, 0, n)]) -= xi * zeta;
        M(i, coord2idx[make_mytuple(u + q, v, n)]) -= xi * zeta;
        M(i, coord2idx[make_mytuple(0, v + u, n)]) -= eta * zeta;
        M(i, coord2idx[make_mytuple(u, v + q, n)]) -= eta * zeta;
        M(i, coord2idx[make_mytuple(0, v, n)]) -= rho * zeta;
        M(i, coord2idx[make_mytuple(u, 0, n)]) -= rho * zeta;
    
        M(i, coord2idx[make_mytuple(n, 0, 0)]) += xi * (1 - zeta);
        M(i, coord2idx[make_mytuple(0, n, 0)]) += eta * (1 - zeta);
        M(i, coord2idx[make_mytuple(0, 0, 0)]) += rho * (1 - zeta);
        M(i, coord2idx[make_mytuple(n, 0, n)]) += xi * zeta;
        M(i, coord2idx[make_mytuple(0, n, n)]) += eta * zeta;
        M(i, coord2idx[make_mytuple(0, 0, n)]) += rho * zeta;
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementPyramid(int order)
    {
      if(order < 3) return fullMatrix<double>(0, 0);
    
      const int szInt = (order - 2) * ((order - 2) + 1) * (2 * (order - 2) + 1) / 6;
      const int szComp =
        (order + 1) * ((order + 1) + 1) * (2 * (order + 1) + 1) / 6;
      const int szInc = szComp - szInt;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsPyramid(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<mytuple, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        int w = order - coordinates(i, 2);
        coord2idx[make_mytuple(u, v, w)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        int w = order - coordinates(szInc + i, 2);
        int q = n - u - w;
        int r = n - v - w;
        double xi = (double)u / n;
        double eta = (double)v / n;
        double rho = (double)q / n;
        double tau = (double)r / n;
        double xip = (double)u / (n - w);
        double etap = (double)v / (n - w);
    
        // lateral faces
        M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - xip;
        M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - etap;
        M(i, coord2idx[make_mytuple(n - w, v, w)]) += xip;
        M(i, coord2idx[make_mytuple(u, n - w, w)]) += etap;
        // basis face
        M(i, coord2idx[make_mytuple(u, v, 0)]) += rho * tau;
        M(i, coord2idx[make_mytuple(u + w, v, 0)]) += xi * tau;
        M(i, coord2idx[make_mytuple(u + w, v + w, 0)]) += xi * eta;
        M(i, coord2idx[make_mytuple(u, v + w, 0)]) += rho * eta;
        // vertical edges
        M(i, coord2idx[make_mytuple(0, 0, w)]) -= (1 - xip) * (1 - etap);
        M(i, coord2idx[make_mytuple(n - w, 0, w)]) -= xip * (1 - etap);
        M(i, coord2idx[make_mytuple(n - w, n - w, w)]) -= xip * etap;
        M(i, coord2idx[make_mytuple(0, n - w, w)]) -= (1 - xip) * etap;
        // basis edges
        M(i, coord2idx[make_mytuple(0, v, 0)]) -= rho * tau;
        M(i, coord2idx[make_mytuple(u, 0, 0)]) -= rho * tau;
        M(i, coord2idx[make_mytuple(u + w, 0, 0)]) -= xi * tau;
        M(i, coord2idx[make_mytuple(n, v, 0)]) -= xi * tau;
        M(i, coord2idx[make_mytuple(n, v + w, 0)]) -= xi * eta;
        M(i, coord2idx[make_mytuple(u + w, n, 0)]) -= xi * eta;
        M(i, coord2idx[make_mytuple(u, n, 0)]) -= rho * eta;
        M(i, coord2idx[make_mytuple(0, v + w, 0)]) -= rho * eta;
        // basis corners
        M(i, coord2idx[make_mytuple(0, 0, 0)]) += rho * tau;
        M(i, coord2idx[make_mytuple(n, 0, 0)]) += xi * tau;
        M(i, coord2idx[make_mytuple(n, n, 0)]) += xi * eta;
        M(i, coord2idx[make_mytuple(0, n, 0)]) += rho * eta;
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementTriangleLinear(int order,
                                                                      int dir)
    {
      if(order < 3) return fullMatrix<double>(0, 0);
    
      const int szInc = 3 * order;
      const int szComp = (order + 1) * (order + 2) / 2;
      const int szInt = szComp - szInc;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsTriangle(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<std::pair<int, int>, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        coord2idx[std::make_pair(u, v)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        double mu;
    
        switch(dir) {
        case 0:
          mu = (double)u / (n - v);
          M(i, coord2idx[std::make_pair(0, v)]) += 1 - mu;
          M(i, coord2idx[std::make_pair(n - v, v)]) += mu;
          break;
        case 1:
          mu = (double)v / (u + v);
          M(i, coord2idx[std::make_pair(u + v, 0)]) += 1 - mu;
          M(i, coord2idx[std::make_pair(0, u + v)]) += mu;
          break;
        default:
        case 2:
          mu = (double)v / (n - u);
          M(i, coord2idx[std::make_pair(u, 0)]) += 1 - mu;
          M(i, coord2idx[std::make_pair(u, n - u)]) += mu;
          break;
        }
      }
    
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementQuadrangleLinear(int order)
    {
      if(order < 2) return fullMatrix<double>(0, 0);
    
      const int szInc = 4 * order;
      const int szComp = (order + 1) * (order + 1);
      const int szInt = szComp - szInc;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsQuadrangle(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<std::pair<int, int>, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        coord2idx[std::make_pair(u, v)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        double eta = (double)v / n;
    
        M(i, coord2idx[std::make_pair(u, 0)]) += 1 - eta;
        M(i, coord2idx[std::make_pair(u, n)]) += eta;
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementTetrahedronLinear(int order,
                                                                         int dir)
    {
      // 'dir' corresponds to the number of edge
      if(order < 4) return fullMatrix<double>(0, 0);
    
      const int szInt = (order - 3) * (order - 2) * (order - 1) / 6;
      const int szComp = (order + 1) * (order + 2) * (order + 3) / 6;
      const int szInc = szComp - szInt;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsTetrahedron(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<mytuple, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        int w = coordinates(i, 2);
        coord2idx[make_mytuple(u, v, w)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        int w = coordinates(szInc + i, 2);
        double mu;
    
        switch(dir) {
        case 0:
          mu = (double)u / (n - v - w);
          M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(n - v - w, v, w)]) += mu;
          break;
        case 2:
          mu = (double)v / (n - u - w);
          M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(u, n - u - w, w)]) += mu;
          break;
        default:
        case 3:
          mu = (double)w / (n - u - v);
          M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(u, v, n - u - v)]) += mu;
          break;
        case 1:
          mu = (double)v / (u + v);
          M(i, coord2idx[make_mytuple(u + v, 0, w)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(0, u + v, w)]) += mu;
          break;
        case 4:
          mu = (double)w / (v + w);
          M(i, coord2idx[make_mytuple(u, v + w, 0)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(u, 0, v + w)]) += mu;
          break;
        case 5:
          mu = (double)w / (u + w);
          M(i, coord2idx[make_mytuple(u + w, v, 0)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(0, v, u + w)]) += mu;
          break;
        }
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementHexahedronLinear(int order,
                                                                        int dir)
    {
      if(order < 2) return fullMatrix<double>(0, 0);
    
      const int szInt = (order - 1) * (order - 1) * (order - 1);
      const int szComp = (order + 1) * (order + 1) * (order + 1);
      const int szInc = szComp - szInt;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsHexahedron(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<mytuple, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        int w = coordinates(i, 2);
        coord2idx[make_mytuple(u, v, w)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        int w = coordinates(szInc + i, 2);
        double eta;
    
        switch(dir) {
        case 0:
          eta = (double)u / n;
          M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - eta;
          M(i, coord2idx[make_mytuple(n, v, w)]) += eta;
          break;
        case 1:
          eta = (double)v / n;
          M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - eta;
          M(i, coord2idx[make_mytuple(u, n, w)]) += eta;
          break;
        default:
        case 2:
          eta = (double)w / n;
          M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - eta;
          M(i, coord2idx[make_mytuple(u, v, n)]) += eta;
          break;
        }
      }
      return M;
    }
    
    fullMatrix<double> gmshGenerateInnerVertexPlacementPrismLinear(int order,
                                                                   int dir)
    {
      if(order < 3) return fullMatrix<double>(0, 0);
    
      const int szInt = (order - 1) * (order - 2) * (order - 1) / 2;
      const int szComp = (order + 1) * (order + 1) * (order + 2) / 2;
      const int szInc = szComp - szInt;
    
      fullMatrix<double> monomials = gmshGenerateMonomialsPrism(order, false);
      fullMatrix<int> coordinates;
      double2int(monomials, coordinates);
    
      std::map<mytuple, int> coord2idx;
      for(int i = 0; i < szInc; ++i) {
        int u = coordinates(i, 0);
        int v = coordinates(i, 1);
        int w = coordinates(i, 2);
        coord2idx[make_mytuple(u, v, w)] = i;
      }
    
      int &n = order;
      fullMatrix<double> M(szInt, szInc, true);
      for(int i = 0; i < szInt; ++i) {
        int u = coordinates(szInc + i, 0);
        int v = coordinates(szInc + i, 1);
        int w = coordinates(szInc + i, 2);
        double eta, mu;
    
        switch(dir) {
        case 0:
          mu = (double)u / (n - v);
          M(i, coord2idx[make_mytuple(0, v, w)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(n - v, v, w)]) += mu;
          break;
        case 1:
          mu = (double)v / (u + v);
          M(i, coord2idx[make_mytuple(u + v, 0, w)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(0, u + v, w)]) += mu;
          break;
        case 2:
          mu = (double)v / (n - u);
          M(i, coord2idx[make_mytuple(u, 0, w)]) += 1 - mu;
          M(i, coord2idx[make_mytuple(u, n - u, w)]) += mu;
          break;
        default:
        case 3:
          eta = (double)w / n;
          M(i, coord2idx[make_mytuple(u, v, 0)]) += 1 - eta;
          M(i, coord2idx[make_mytuple(u, v, n)]) += eta;
          break;
        }
      }
      return M;
    }
    
    #undef make_mytuple