Skip to content
Snippets Groups Projects
Select Git revision
  • 72b930c272dca44e83645585075e788a3b180370
  • master default
  • cgnsUnstructured
  • partitioning
  • poppler
  • HighOrderBLCurving
  • gmsh_3_0_4
  • gmsh_3_0_3
  • gmsh_3_0_2
  • gmsh_3_0_1
  • gmsh_3_0_0
  • gmsh_2_16_0
  • gmsh_2_15_0
  • gmsh_2_14_1
  • gmsh_2_14_0
  • gmsh_2_13_2
  • gmsh_2_13_1
  • gmsh_2_12_0
  • gmsh_2_11_0
  • gmsh_2_10_1
  • gmsh_2_10_0
  • gmsh_2_9_3
  • gmsh_2_9_2
  • gmsh_2_9_1
  • gmsh_2_9_0
  • gmsh_2_8_6
26 results

ReferenceSpace.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    ReferenceSpace.cpp 31.80 KiB
    #include <algorithm>
    #include <fstream>
    #include <sstream>
    #include <cstring>
    #include <set>
    
    #include "Exception.h"
    #include "Numeric.h"
    #include "ReferenceSpace.h"
    
    using namespace std;
    
    ReferenceSpace::ReferenceSpace(void){
      // Defining Ref Edge and Face in //
      // Derived Class                 //
    }
    
    ReferenceSpace::ReferenceSpace(const std::string& path){
      init(path);
    }
    
    ReferenceSpace::~ReferenceSpace(void){
      delete pTree;
    }
    
    void ReferenceSpace::init(void){
      // Tree //
      vector<size_t> vertexSeq(nVertex);
    
      for(size_t i = 0; i < nVertex; i++)
        vertexSeq[i] = i;
    
      pTree = new PermutationTree(vertexSeq);
    
      // Cyclic Permutation //
      list<size_t>          listOfTrueReferenceSpace;
      list<vector<size_t> > listOfRefNodeIndexPermutation;
      list<vector<size_t> > listOfReverseNodeIndexPermutation;
    
      findCyclicPermutation(listOfTrueReferenceSpace,
                            listOfRefNodeIndexPermutation,
                            listOfReverseNodeIndexPermutation);
    
      // Iterators //
      list<size_t>::iterator refSpaceIt = listOfTrueReferenceSpace.begin();
    
      // Reference Spaces Node Id //
      const size_t nRefSpace = listOfTrueReferenceSpace.size();
      refSpaceNodeId.resize(nRefSpace);
    
      for(size_t i = 0; i < nRefSpace; i++, refSpaceIt++){
        refSpaceNodeId[i].resize(nVertex);
    
        pTree->fillWithPermutation(*refSpaceIt, refSpaceNodeId[i]);
      }
    
      // Compress Tag, so that they refer to refSpaceNodeId indices
      // insted of leafId of pTree.
      pTree->compressTag();
    
      // ABC space to UVW space shape function index mapping //
      const size_t nPermutation = pTree->getNPermutation();
      ABCtoUVWIndex.resize(nPermutation);
      ABCtoUVWIndex.assign(listOfReverseNodeIndexPermutation.begin(),
                           listOfReverseNodeIndexPermutation.end());
    
      // UVW space to ABC space shape function index mapping //
      UVWtoABCIndex.resize(nPermutation);
      UVWtoABCIndex.assign(listOfRefNodeIndexPermutation.begin(),
                           listOfRefNodeIndexPermutation.end());
    
      // Ordered Edges & Faces Node Index//
      getOrderedEdge();
      getOrderedFace();
    }
    
    void ReferenceSpace::init(const char* stream){
      // nVector
      size_t offset = 0;
    
      memcpy(&nVertex, stream + offset, sizeof(size_t));
      offset += sizeof(size_t);
    
      // Fef {Edge, Face} Node Idx
      offset += unserialize(stream + offset, refEdgeNodeIdx);
      offset += unserialize(stream + offset, refFaceNodeIdx);
    
      // Permutation Tree
      size_t tSize; // Tree Size
    
      memcpy(&tSize, stream + offset, sizeof(size_t));
      offset += sizeof(size_t);
    
      pTree = new PermutationTree(stream + offset);
      offset += tSize;
    
      // RefSpace Node Id
      offset += unserialize(stream + offset, refSpaceNodeId);
    
      // Index mapping
      offset += unserialize(stream + offset, ABCtoUVWIndex);
      offset += unserialize(stream + offset, UVWtoABCIndex);
    
      // Oredered {Edge, Face} Node Idx
      offset += unserialize(stream + offset, orderedEdgeNodeIdx);
      offset += unserialize(stream + offset, orderedFaceNodeIdx);
    }
    
    void ReferenceSpace::init(const std::string& path){
      // Read file //
      // Open Stream
      ifstream input;
      input.exceptions(std::ifstream::failbit | std::ifstream::badbit);
      input.open(path.c_str(), std::ifstream::binary);
    
      // Get size of stream (go to stream end)
      input.seekg(0, std::ifstream::end);
      const size_t size = input.tellg();
    
      // Reset stream possition
      input.seekg(0, std::ifstream::beg);
    
      // Alloc byte stream & Read file
      char* stream = new char[size];
      input.read(stream, size);
    
      // Init from stream
      init(stream);
    
      // Free stream
      delete[] stream;
    }
    
    void ReferenceSpace::
    findCyclicPermutation(list<size_t>&          listOfTrueReferenceSpace,
                          list<vector<size_t> >& listOfRefNodeIndexPermutation,
                          list<vector<size_t> >& listOfReverseNodeIndexPermutation){
      // Alloc Some Data //
      const size_t nPerm = pTree->getNPermutation();
    
      vector<size_t> pTest(nVertex);
      vector<size_t> pRef(nVertex);
    
      list<size_t>::iterator it;
      list<size_t>::iterator end;
      triplet match;
    
      // Values For unPermutedIndex //
      vector<size_t> unPermutedIndex(nVertex);
    
      for(size_t i = 0; i < nVertex; i++)
        unPermutedIndex[i] = i;
    
      // Find Cyclic Permutation //
      for(size_t i = 0; i < nPerm; i++){
        // No match
        match.first = false;
    
        // Get Permutation 'i'
        pTree->fillWithPermutation(i, pTest);
    
        // Test it with already found Reference Space
        it  = listOfTrueReferenceSpace.begin();
        end = listOfTrueReferenceSpace.end();
    
        while(it != end && !match.first){
          // Take Reference Space 'it'
          pTree->fillWithPermutation(*it, pRef);
    
          // Look if it matches
          match = isCyclicPermutation(pTest, pRef);
    
          // If not, go to next Reference Space
          if(!match.first)
            it++;
        }
    
        // If no Reference Space is found
        //   --> this Permutation is a new Reference Space
        if(!match.first){
          listOfTrueReferenceSpace.push_back(i);
    
          listOfRefNodeIndexPermutation.push_back(unPermutedIndex);
          listOfReverseNodeIndexPermutation.push_back(unPermutedIndex);
    
          pTree->addTagToPermutation(i, i);
        }
    
        // If a ReferenceSpace is found, and the index permutations
        else{
          listOfRefNodeIndexPermutation.push_back(match.second);
          listOfReverseNodeIndexPermutation.push_back(match.third);
    
          pTree->addTagToPermutation(i, *it);
        }
      }
    }
    
    ReferenceSpace::triplet ReferenceSpace::
    isCyclicPermutation(const vector<size_t>& pTest,
                        const vector<size_t>& pRef) const{
    
      // Triplet to return
      triplet tri;
      /*
      // If no Face, we have no Cyclic Permutation
      if(!refFaceNodeIdx.size()){
        tri.first  = false;
        tri.second = vector<size_t>(0);
        tri.third  = vector<size_t>(0);
    
        return tri;
      }
    
      // Node IDs of Reference Space first Face
      const size_t   nNodeInFaceZero = refFaceNodeIdx[0].size();
      vector<size_t> refNode(nNodeInFaceZero);
    
      for(size_t i = 0; i < nNodeInFaceZero; i++)
        refNode[i] = pRef[refFaceNodeIdx[0][i]];
    
      // Corresponding Face in Test Permutation
      size_t testFaceId = findCorrespondingFace(refNode, pTest);
    
      // If no corresponding face found, we have no Cyclic Permutation
      if(testFaceId == (size_t)(-1)){
        tri.first  = false;
        tri.second = vector<size_t>(0);
        tri.third  = vector<size_t>(0);
    
        return tri;
      }
    
      // Node IDs of Test Permutation correspnding Face
      const size_t   nNodeInFaceTest = refFaceNodeIdx[testFaceId].size();
      vector<size_t> testNode(nNodeInFaceTest);
    
      for(size_t i = 0; i < nNodeInFaceTest; i++)
        testNode[i] = pTest[refFaceNodeIdx[testFaceId][i]];
    
      // Test if we have a face cyclic permutation
      size_t isCyclic = isFacePermutation(refNode, testNode);
      */
      // Test if we have the same connectivity
      bool isSameConnectivity = isSameEdge(pTest, pRef);
    
      if(isSameConnectivity){
        tri.first  = true;
        tri.second = getRefIndexPermutation(pRef, pTest);
        tri.third  = getReverseIndexPermutation(pRef, pTest);
      }
    
      else{
        tri.first  = false;
        tri.second = vector<size_t>(0);
        tri.third  = vector<size_t>(0);
      }
    
      return tri;
    }
    /*
    size_t ReferenceSpace::findCorrespondingFace(const vector<size_t>& face,
                                                 const vector<size_t>& node) const{
      // Init Stuff //
      const size_t nFace    = refFaceNodeIdx.size();
      const size_t faceSize = face.size();
      bool match = false;
      size_t f = 0;
    
      vector<size_t> testFace(faceSize);
      size_t         nNodeInFace;
    
      // Test All Face until match
      while(!match && f < nFace){
        nNodeInFace = refFaceNodeIdx[f].size();
    
        if(nNodeInFace == faceSize){
          // Get face f nodes
          for(size_t i = 0; i < faceSize; i++)
            testFace[i] = node[refFaceNodeIdx[f][i]];
    
          // Look if match
          match = haveSameNode(testFace, face);
        }
    
        if(!match)
          f++;
      }
    
      // If a valid face is found return it
      if(f < nFace)
        return f;
    
      // Else return (size_t)(-1)
      else
        return -1;
    }
    
    bool ReferenceSpace::haveSameNode(const vector<size_t>& face0,
                                      const vector<size_t>& face1){
      // Check Sizes
      const size_t size = face0.size();
    
      if(size != face1.size())
        return false;
    
      // Copy vector
      vector<size_t> cFace0(face0);
      vector<size_t> cFace1(face1);
    
      // Sort both copy vectors
      sort(cFace0.begin(), cFace0.end());
      sort(cFace1.begin(), cFace1.end());
    
      // Check if elements are the same and at same position
      bool match = true;
    
      for(size_t i = 0; i < size && match; i++)
        match = (cFace0[i] == cFace1[i]);
    
      return match;
    }
    
    bool ReferenceSpace::isFacePermutation(const vector<size_t>& refNode,
                                           const vector<size_t>& testNode){
      // Check Size
      const size_t size = refNode.size();
    
      if(size != testNode.size())
        return false;
    
      // Find refNode[0] in testNode
      size_t offset = 0;
    
      while(refNode[0] != testNode[offset] && offset < size)
        offset++;
    
      // If refNode[0] not found in testNode, there is no face permutation
      if(offset == size)
        return false;
    
      // Check refNode[i] and testNode[(i + offset) % size]
      bool match = true;
    
      for(size_t i = 0; i < size && match; i++)
        match = (refNode[i] == testNode[(i + offset) % size]);
    
      return match;
    }
    */
    bool ReferenceSpace::isSameEdge(const std::vector<size_t>& pTest,
                                    const std::vector<size_t>& pRef) const{
      // Set of Reference Edges
      const size_t nEdge = refEdgeNodeIdx.size();
      std::set<std::vector<size_t>, EdgeComparator> refEdge;
      std::vector<size_t> tmp(2);
    
      for(size_t e = 0; e < nEdge; e++){
        tmp[0] = pRef[refEdgeNodeIdx[e][0]];
        tmp[1] = pRef[refEdgeNodeIdx[e][1]];
        refEdge.insert(tmp);
      }
    
      // Compare Test Edges
      bool match = true;
      for(size_t e = 0; e < nEdge && match; e++){
        tmp[0] = pTest[refEdgeNodeIdx[e][0]];
        tmp[1] = pTest[refEdgeNodeIdx[e][1]];
    
        match = (refEdge.find(tmp) != refEdge.end());
      }
    
      return match;
    }
    
    bool ReferenceSpace::EdgeComparator::operator()(const std::vector<size_t>& a,
                                                    const std::vector<size_t>& b){
      size_t maxA, minA;
      size_t maxB, minB;
    
      if(a[0] < a[1]){
        maxA = a[1];
        minA = a[0];
      }
    
      else{
        maxA = a[0];
        minA = a[1];
      }
    
      if(b[0] < b[1]){
        maxB = b[1];
        minB = b[0];
      }
    
      else{
        maxB = b[0];
        minB = b[1];
      }
    
      return
        ((minA != minB) && (minA < minB)) ||
        ((minA == minB) && (maxA < maxB));
    }
    
    vector<size_t>
    ReferenceSpace::getRefIndexPermutation(const vector<size_t>& ref,
                                           const vector<size_t>& test) const{
      const size_t size = ref.size();
      vector<size_t> idxVec(ref.size());
      size_t idx;
    
      for(size_t i = 0; i < size; i++){
        idx = 0;
    
        while(test[i] != ref[idx])
          idx++;
    
        idxVec[i] = idx;
      }
    
      return idxVec;
    }
    
    vector<size_t>
    ReferenceSpace::getReverseIndexPermutation(const vector<size_t>& ref,
                                               const vector<size_t>& test) const{
      const size_t size = ref.size();
      vector<size_t> idxVec(ref.size());
      size_t idx;
    
      for(size_t i = 0; i < size; i++){
        idx = 0;
    
        while(test[idx] != ref[i])
          idx++;
    
        idxVec[i] = idx;
      }
    
      return idxVec;
    }
    
    void ReferenceSpace::getOrderedEdge(void){
      // Fill orderedEdgeNodeIdx[s][e]                         //
      // (for all reference space 's' and edge 'e') such that: //
      //                                                       //
      // refSpaceNodeId[s][orderedEdgeNodeIdx[s][e][0]] <      //
      // refSpaceNodeId[s][orderedEdgeNodeIdx[s][e][1]]        //
    
      // Alloc
      const size_t nEdge     = refEdgeNodeIdx.size();
      const size_t nRefSpace = refSpaceNodeId.size();
    
      orderedEdgeNodeIdx.resize(nRefSpace);
    
      // Populate Edge
      for(size_t s = 0; s < nRefSpace; s++){
        orderedEdgeNodeIdx[s].resize(nEdge);
    
        for(size_t e = 0; e < nEdge; e++){
          orderedEdgeNodeIdx[s][e].resize(2);
    
          orderRefEntityForGivenRefSpace(refEdgeNodeIdx[e],
                                         refSpaceNodeId[s],
                                         orderedEdgeNodeIdx[s][e]);
        }
      }
    }
    
    void ReferenceSpace::getOrderedFace(void){
      // Fill orderedEdgeNodeIdx[s][f]                                //
      // (for all reference space 's' and face 'f') such that:        //
      //                                                              //
      // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][0]] <             //
      // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][1]] <             //
      // refSpaceNodeId[s][orderedFaceNodeIdx[s][f][2]] <             //
      // (refSpaceNodeId[s][orderedFaceNodeIdx[s][f][3]]) (Quad Face) //
      //                                                              //
      // If we have a Quad face, a correction is needed such that:    //
      //                                                              //
      // orderedFaceNodeIdx[s][f][2]                                  //
      // is *opposite* to                                             //
      // orderedFaceNodeIdx[s][f][0]                                  //
    
      // Alloc
      const size_t nFace     = refFaceNodeIdx.size();
      const size_t nRefSpace = refSpaceNodeId.size();
    
      orderedFaceNodeIdx.resize(nRefSpace);
    
      // Populate Face
      for(size_t s = 0; s < nRefSpace; s++){
        orderedFaceNodeIdx[s].resize(nFace);
    
        for(size_t f = 0; f < nFace; f++){
          size_t nNodeInFace = refFaceNodeIdx[f].size();
    
          orderedFaceNodeIdx[s][f].resize(nNodeInFace);
    
          orderRefEntityForGivenRefSpace(refFaceNodeIdx[f],
                                         refSpaceNodeId[s],
                                         orderedFaceNodeIdx[s][f]);
          if(nNodeInFace == 4)
            correctQuadFaceNodeIdx(f, orderedFaceNodeIdx[s][f]);
        }
      }
    }
    
    void ReferenceSpace::
    orderRefEntityForGivenRefSpace(vector<size_t>& refEntityNodeIdx,
                                   vector<size_t>& refSpaceNodeId,
                                   vector<size_t>& orderedEntityNodeIdx){
      // Get Size
      const size_t size = orderedEntityNodeIdx.size();
    
      // Ref Entity Node *ID*
      vector<pair<size_t, size_t> > refEntityNodeId(size);
    
      for(size_t i = 0; i < size; i++){
        refEntityNodeId[i].first  = i;                             // Node Id INDEX
        refEntityNodeId[i].second = refSpaceNodeId[refEntityNodeIdx[i]]; // Node Id
      }
    
      // Sort it with repsect to second entry (Node ID)
      std::sort(refEntityNodeId.begin(), refEntityNodeId.end(), sortPredicate);
    
      // Populate orderedEntityNodeIdx (Usign sorted Node ID old Index)
      for(size_t i = 0; i < size; i++)
        orderedEntityNodeIdx[i] = refEntityNodeIdx[refEntityNodeId[i].first];
    }
    
    void ReferenceSpace::
    correctQuadFaceNodeIdx(size_t faceId,
                           vector<size_t>& correctedQuadFaceNodeIdx){
      // Get :                       //
      // correctedQuadFaceNodeIdx[2] //
      // is *opposite* to            //
      // correctedQuadFaceNodeIdx[0] //
    
      // Get opposite Node //
      // Find node correctedQuadFaceNodeIdx[0] in refFaceNodeId[faceId]
      size_t refIdx = 0;
    
      while((refIdx < 4)
            &&
            (refFaceNodeIdx[faceId][refIdx] != correctedQuadFaceNodeIdx[0])){
    
        refIdx++;
      }
    
      if(refIdx == 4)
        throw Exception("Error in correctQuadFaceNodeIdx: vertex 0 not found");
    
      // Get opposite Node
      size_t opposite = refFaceNodeIdx[faceId][(refIdx + 2) % 4];
    
      // If correction needed, do it //
      if(correctedQuadFaceNodeIdx[2] != opposite){
        // Find opposite node index
        size_t tmp;
        size_t idx = 1;
    
        while(correctedQuadFaceNodeIdx[idx] != opposite)
          idx++;
    
        // Swap correctedQuadFaceNodeIdx[2] and correctedQuadFaceNodeIdx[idx]
        tmp = correctedQuadFaceNodeIdx[2];
        correctedQuadFaceNodeIdx[2]   = opposite;
        correctedQuadFaceNodeIdx[idx] = tmp;
      }
    }
    
    size_t ReferenceSpace::getPermutationIdx(const MElement& element) const{
      // Get Primary Vertices Global ID //
      vector<pair<size_t, size_t> > vertexGlobalId(nVertex);
    
      for(size_t i = 0; i < nVertex; i++){
        vertexGlobalId[i].first  = i;
        vertexGlobalId[i].second = element.getVertex(i)->getNum();
      }
    
      // Sort Them with repsect to Vertex Global ID //
      //                 (vertex[i].second->getNum) //
      std::sort(vertexGlobalId.begin(), vertexGlobalId.end(), sortPredicate);
    
      // Reduce Vertex Global ID //
      vector<size_t> vertexReducedId(nVertex);
    
      for(size_t i = 0; i < nVertex; i++)
        vertexReducedId[vertexGlobalId[i].first] = i;
    
      // Tree Lookup //
      try{
        return pTree->getPermutationId(vertexReducedId);
      }
    
      catch(...){
        throw Exception("Cannot Find Reference Space for Element %d",
                        element.getNum());
      }
    }
    
    void ReferenceSpace::mapFromABCtoUVW(const MElement& element,
                                         double a, double b, double c,
                                         double uvw[3]) const{
      // Get Index Permutation
      const size_t permutationIdx = getPermutationIdx(element);
    
      // UVW node coordinate in ABC
      double** uvwNode = new double*[nVertex];
    
      for(size_t i = 0; i < nVertex; i++)
        uvwNode[i] = new double[3];
    
      for(size_t i = 0; i < nVertex; i++)
        element.getNode(i,
                        uvwNode[UVWtoABCIndex[permutationIdx][i]][0],
                        uvwNode[UVWtoABCIndex[permutationIdx][i]][1],
                        uvwNode[UVWtoABCIndex[permutationIdx][i]][2]);
    
      // ABC (order 1) grad shape functions
      double* phiABC = new double[nVertex];
      element.getShapeFunctions(a, b, c, phiABC, 1);
    
      // Map From ABC to UVW //
      uvw[0] = 0;
      for(size_t i = 0; i < nVertex; i++)
        uvw[0] += uvwNode[i][0] * phiABC[i];
    
      uvw[1] = 0;
      for(size_t i = 0; i < nVertex; i++)
        uvw[1] += uvwNode[i][1] * phiABC[i];
    
      uvw[2] = 0;
      for(size_t i = 0; i < nVertex; i++)
        uvw[2] += uvwNode[i][2] * phiABC[i];
    
      // Free //
      delete[] phiABC;
    
      for(size_t i = 0; i < nVertex; i++)
        delete[] uvwNode[i];
    
      delete[] uvwNode;
    }
    
    void ReferenceSpace::mapFromABCtoXYZ(const MElement& element,
                                         double a, double b, double c,
                                         double xyz[3]) const{
      // Map From ABC to UVW
      double uvw[3];
      mapFromABCtoUVW(element, a, b, c, uvw);
    
      // Map From UVW to XYZ
      SPoint3 pxyz;
      element.pnt(uvw[0], uvw[1], uvw[2], pxyz);
    
      xyz[0] = pxyz.x();
      xyz[1] = pxyz.y();
      xyz[2] = pxyz.z();
    }
    
    void ReferenceSpace::mapFromUVWtoABC(const MElement& element,
                                         double u, double v, double w,
                                         double abc[3]) const{
      // Element Permutation Index
      const size_t permutationIdx = getPermutationIdx(element);
    
      // ABC node coordinate in UVW
      double** abcNode = new double*[nVertex];
    
      for(size_t i = 0; i < nVertex; i++)
        abcNode[i] = new double[3];
    
      for(size_t i = 0; i < nVertex; i++)
        element.getNode(i,
                        abcNode[ABCtoUVWIndex[permutationIdx][i]][0],
                        abcNode[ABCtoUVWIndex[permutationIdx][i]][1],
                        abcNode[ABCtoUVWIndex[permutationIdx][i]][2]);
    
      // UVW (order 1) shape functions
      double* phiUVW = new double[nVertex];
      element.getShapeFunctions(u, v, w, phiUVW, 1);
    
      // Map From UVW to ABC
      abc[0] = 0;
      for(size_t i = 0; i < nVertex; i++)
        abc[0] += abcNode[i][0] * phiUVW[i];
    
      abc[1] = 0;
      for(size_t i = 0; i < nVertex; i++)
        abc[1] += abcNode[i][1] * phiUVW[i];
    
      abc[2] = 0;
      for(size_t i = 0; i < nVertex; i++)
        abc[2] += abcNode[i][2] * phiUVW[i];
    
      // Free
      delete[] phiUVW;
    
      for(size_t i = 0; i < nVertex; i++)
        delete[] abcNode[i];
    
      delete[] abcNode;
    }
    
    void ReferenceSpace::mapFromXYZtoABC(const MElement& element,
                                         double x, double y, double z,
                                         double abc[3]) const{
      // Get UVW coordinate //
      double xyz[3] = {x, y, z};
      double uvw[3];
    
      element.xyz2uvw(xyz, uvw);
    
      // Get ABC coordinate //
      mapFromUVWtoABC(element, uvw[0], uvw[1], uvw[2], abc);
    }
    
    double ReferenceSpace::getJacobian(const MElement& element,
                                       double a, double b, double c,
                                       fullMatrix<double>& jac) const{
      // ABC to UVW Jacobian //
      // Get Index Permutation
      const size_t permutationIdx = getPermutationIdx(element);
    
      // UVW node coordinate
      double** uvwNode = new double*[nVertex];
    
      for(size_t i = 0; i < nVertex; i++)
        uvwNode[i] = new double[3];
    
      for(size_t i = 0; i < nVertex; i++)
        element.getNode(i,
                        uvwNode[UVWtoABCIndex[permutationIdx][i]][0],
                        uvwNode[UVWtoABCIndex[permutationIdx][i]][1],
                        uvwNode[UVWtoABCIndex[permutationIdx][i]][2]);
    
      // ABC (order 1) grad shape functions
      double phiABC[1256][3]; // Cannot be dynamicaly allocated since
                              // GMSH wants a double[][3] for phiABC !!!
    
      element.getGradShapeFunctions(a, b, c, phiABC, 1);
    
      // Jacobian
      fullMatrix<double> jacABCtoUVW(3, 3);
      jacABCtoUVW.setAll(0);
    
      for(size_t i = 0; i < nVertex; i++){
        jacABCtoUVW(0, 0) += uvwNode[i][0] * phiABC[i][0];
        jacABCtoUVW(0, 1) += uvwNode[i][1] * phiABC[i][0];
        jacABCtoUVW(0, 2) += uvwNode[i][2] * phiABC[i][0];
      }
    
      for(size_t i = 0; i < nVertex; i++){
        jacABCtoUVW(1, 0) += uvwNode[i][0] * phiABC[i][1];
        jacABCtoUVW(1, 1) += uvwNode[i][1] * phiABC[i][1];
        jacABCtoUVW(1, 2) += uvwNode[i][2] * phiABC[i][1];
      }
    
      for(size_t i = 0; i < nVertex; i++){
        jacABCtoUVW(2, 0) += uvwNode[i][0] * phiABC[i][2];
        jacABCtoUVW(2, 1) += uvwNode[i][1] * phiABC[i][2];
        jacABCtoUVW(2, 2) += uvwNode[i][2] * phiABC[i][2];
      }
    
      // Regularize Jacobian
      regularize(element.getDim(), jacABCtoUVW);
    
      // Free
      for(size_t i = 0; i < nVertex; i++)
        delete[] uvwNode[i];
    
      delete[] uvwNode;
    
      // Map ABC coordinate to UVW point //
      double uvw[3];
      mapFromABCtoUVW(element, a, b, c, uvw);
    
      // UVW to XYZ Jacobian + Determinant                      //
      //  NB: Volume is not modified when we go from ABC to UVW //
      //       --> Determinant unchanged                        //
      fullMatrix<double> jacUVWtoXYZ(3, 3);
      double det = element.getJacobian(uvw[0], uvw[1], uvw[2], jacUVWtoXYZ);
    
      // Product of the two Jacobians & Return //
      // Do a naive gemm, so that we do not encounter nested threads
      // (limitation of OpenBLAS)
      jac.gemm_naive(jacABCtoUVW, jacUVWtoXYZ, 1, 0);
      return det;
    }
    
    void ReferenceSpace::regularize(size_t dim, fullMatrix<double>& jac){
      double a[3];
      double b[3];
      double c[3];
    
      switch(dim){
      case 0:
        jac(0, 0) = jac(1, 1) = jac(2, 2) = 1.0;
        jac(0, 1) = jac(1, 0) = jac(2, 0) = 0.0;
        jac(0, 2) = jac(1, 2) = jac(2, 1) = 0.0;
        break;
    
      case 1:
        a[0] = jac(0, 0);
        a[1] = jac(0, 1);
        a[2] = jac(0, 2);
    
       if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
           (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))){
          b[0] =  a[1];
          b[1] = -a[0];
          b[2] =  0.;
        }
    
        else{
          b[0] =  0.;
          b[1] =  a[2];
          b[2] = -a[1];
        }
    
        norme(b);
        prodve(a, b, c);
        norme(c);
    
        jac(1, 0) = b[0]; jac(1, 1) = b[1]; jac(1, 2) = b[2];
        jac(2, 0) = c[0]; jac(2, 1) = c[1]; jac(2, 2) = c[2];
        break;
    
      case 2:
        a[0] = jac(0, 0);
        a[1] = jac(0, 1);
        a[2] = jac(0, 2);
    
        b[0] = jac(1, 0);
        b[1] = jac(1, 1);
        b[2] = jac(1, 2);
    
        prodve(a, b, c);
        norme(c);
    
        jac(2, 0) = c[0]; jac(2, 1) = c[1]; jac(2, 2) = c[2];
        break;
    
      case 3:
        break;
      }
    }
    
    string ReferenceSpace::toString(void) const{
      const size_t nRefSpace = refSpaceNodeId.size();
      const size_t nEdge     = refEdgeNodeIdx.size();
      const size_t nFace     = refFaceNodeIdx.size();
      stringstream stream;
    
      // Reference Space Node IDs //
      stream << "Reference Space Node IDs:" << endl;
    
      for(size_t s = 0; s < nRefSpace; s++){
        stream << "  ** Reference Space #" << s << ": [";
    
        for(size_t v = 0; v < nVertex - 1; v++)
          stream << refSpaceNodeId[s][v] << ", ";
    
        stream << refSpaceNodeId[s][nVertex - 1] << "]" << endl;
      }
    
      // Oriented Edges Node Index //
      stream << endl << "Ordered Edge Node Index:" << endl;
    
      for(size_t s = 0; s < nRefSpace; s++){
        stream << "  ** Reference Space #" << s << ":"
               << endl;
    
        for(size_t e = 0; e < nEdge; e++)
          stream << "      -- ["
                 << orderedEdgeNodeIdx[s][e][0] << ", "
                 << orderedEdgeNodeIdx[s][e][1] << "]"
                 << endl;
      }
    
      // Oriented Faces Node Index //
      stream << endl << "Ordered Face Node Index:" << endl;
    
      for(size_t s = 0; s < nRefSpace; s++){
        stream << "  ** Reference Space #" << s << ":"
               << endl;
    
        for(size_t f = 0; f < nFace; f++){
          const size_t nNodeInFace = orderedFaceNodeIdx[s][f].size();
          stream << "      -- [";
    
          for(size_t n = 0; n < nNodeInFace - 1; n++)
            stream << orderedFaceNodeIdx[s][f][n] << ", ";
    
          stream << orderedFaceNodeIdx[s][f][nNodeInFace - 1] << "]"
                 << endl;
        }
      }
    
      return stream.str();
    }
    
    string ReferenceSpace::toLatex(void) const{
      stringstream stream;
    
      stream << "\\documentclass{article}" << endl << endl
             << "\\begin{document}"        << endl
    
             << "\texttt{toLatex} not implemented" << endl
    
             << "\\end{document}"          << endl;
    
      return stream.str();
    }
    
    std::pair<size_t, char*> ReferenceSpace::serialize(void) const{
      // Serialize //
      // refEdgeNodeIdx
      pair<size_t, char*> sRefEdgeIdx = serialize(refEdgeNodeIdx);
    
      // refFaceNodeIdx
      pair<size_t, char*> sRefFaceIdx = serialize(refFaceNodeIdx);
    
      // PermutationTree
      pair<size_t, char*> sPTree = pTree->serialize();
    
      // Reference Space
      pair<size_t, char*> sRSpace = serialize(refSpaceNodeId);
    
      // ABCtoUVW
      pair<size_t, char*> sABCtoUVW = serialize(ABCtoUVWIndex);
    
      // UVWtoUVW
      pair<size_t, char*> sUVWtoABC = serialize(UVWtoABCIndex);
    
      // orderedEdgeNodeIdx
      pair<size_t, char*> sEdgeIdx = serialize(orderedEdgeNodeIdx);
    
      // orderedEdgeNodeIdx
      pair<size_t, char*> sFaceIdx = serialize(orderedFaceNodeIdx);
    
      // Full Serialized Stream //
      // Alloc (+ 1 * sizeof(size_t) for nVertex)
      const size_t size = sizeof(size_t) + sRefEdgeIdx.first
                                         + sRefFaceIdx.first
                                         + sizeof(size_t) // pTree Size
                                         + sPTree.first   // pTree
                                         + sRSpace.first
                                         + sABCtoUVW.first
                                         + sUVWtoABC.first
                                         + sEdgeIdx.first
                                         + sFaceIdx.first;
      char* stream = new char[size];
    
      // Populate
      size_t offset = 0;
    
      memcpy(stream + offset, &nVertex, sizeof(size_t));
      offset += sizeof(size_t);
    
      memcpy(stream + offset, sRefEdgeIdx.second, sRefEdgeIdx.first);
      offset += sRefEdgeIdx.first;
    
      memcpy(stream + offset, sRefFaceIdx.second, sRefFaceIdx.first);
      offset += sRefFaceIdx.first;
    
      memcpy(stream + offset, &sPTree.first, sizeof(size_t)); // pTree size
      offset += sizeof(size_t);
    
      memcpy(stream + offset, sPTree.second, sPTree.first);   // pTree
      offset += sPTree.first;
    
      memcpy(stream + offset, sRSpace.second, sRSpace.first);
      offset += sRSpace.first;
    
      memcpy(stream + offset, sABCtoUVW.second, sABCtoUVW.first);
      offset += sABCtoUVW.first;
    
      memcpy(stream + offset, sUVWtoABC.second, sABCtoUVW.first);
      offset += sABCtoUVW.first;
    
      memcpy(stream + offset, sEdgeIdx.second, sEdgeIdx.first);
      offset += sEdgeIdx.first;
    
      memcpy(stream + offset, sFaceIdx.second, sFaceIdx.first);
      offset += sFaceIdx.first;
    
      // Free Temp
      delete[] sRefEdgeIdx.second;
      delete[] sRefFaceIdx.second;
      delete[] sPTree.second;
      delete[] sRSpace.second;
      delete[] sABCtoUVW.second;
      delete[] sUVWtoABC.second;
      delete[] sEdgeIdx.second;
      delete[] sFaceIdx.second;
    
      // Return
      return pair<size_t, char*>(size, stream);
    }
    
    void ReferenceSpace::serialize(const std::string& path) const{
      // Serialize into byte Stream
      std::pair<size_t, char*> stream  = serialize();
    
      // Write byte stream
      ofstream output;
      output.exceptions(std::ofstream::failbit | std::ofstream::badbit);
      output.open(path.c_str(), std::ofstream::binary);
      output.write(stream.second, stream.first);
      output.close();
    
      // Free
      delete[] stream.second;
    }
    
    std::pair<size_t, char*> ReferenceSpace::
    serialize(const std::vector<std::vector<size_t> >& source){
      // Vector & Stream size //
      size_t nVector   = source.size();
      size_t totalSize = sizeof(size_t); // Init to header size (nVector)
    
      for(size_t i = 0; i < nVector; i++)
        totalSize += (source[i].size() + 1) * sizeof(size_t);
                                    // + 1 for subvector size header
    
      // Alloc Stream //
      char* stream = new char[totalSize];
    
      // Header (nVector) //
      memcpy(stream, &nVector, sizeof(size_t));
    
      // SubVectors //
      size_t tmpSize;
      size_t offset = sizeof(size_t);
    
      for(size_t i = 0; i < nVector; i++){
        // SubVector size
        tmpSize = source[i].size();
        memcpy(stream + offset, &tmpSize, sizeof(size_t));
        offset += sizeof(size_t);
    
        // SubVector stored values
        for(size_t j = 0; j < tmpSize; j++){
          memcpy(stream + offset, &source[i][j], sizeof(size_t));
          offset += sizeof(size_t);
        }
      }
    
      // Return //
      return pair<size_t, char*>(totalSize, stream);
    }
    
    std::pair<size_t, char*> ReferenceSpace::
    serialize(const std::vector<std::vector<std::vector<size_t> > >& source){
      // Vector & Stream size //
      size_t nVector   = source.size();
      size_t totalSize = sizeof(size_t); // Init to header size (nVector)
    
      for(size_t i = 0; i < nVector; i++){
        totalSize += sizeof(size_t); // Subvector size header
    
        for(size_t j = 0; j < source[i].size(); j++)
          totalSize += (source[i][j].size() + 1) * sizeof(size_t);
                                         // + 1 for subsubvector size header
      }
    
      // Alloc Stream //
      char* stream = new char[totalSize];
    
      // Header (nVector) //
      memcpy(stream, &nVector, sizeof(size_t));
    
      // SubVectors //
      size_t tmpSize;
      size_t tmpSubSize;
      size_t offset = sizeof(size_t);
    
      for(size_t i = 0; i < nVector; i++){
        // SubVector size
        tmpSize = source[i].size();
        memcpy(stream + offset, &tmpSize, sizeof(size_t));
        offset += sizeof(size_t);
    
        for(size_t j = 0; j < tmpSize; j++){
          // SubSubVector size
          tmpSubSize = source[i][j].size();
          memcpy(stream + offset, &tmpSubSize, sizeof(size_t));
          offset += sizeof(size_t);
    
          // SubSubVector stored values
          for(size_t k = 0; k < tmpSubSize; k++){
            memcpy(stream + offset, &source[i][j][k], sizeof(size_t));
            offset += sizeof(size_t);
          }
        }
      }
    
      // Return //
      return pair<size_t, char*>(totalSize, stream);
    }
    
    size_t ReferenceSpace::
    unserialize(const char* stream,
                std::vector<std::vector<size_t> >& dest){
      // Vector size //
      size_t nVector;
      size_t offset = 0;
    
      memcpy(&nVector, stream + offset, sizeof(size_t));
      offset += sizeof(size_t);
    
      // Alloc Vector //
      dest.resize(nVector);
    
      // Populate //
      size_t tmpSize;
      size_t tmpVal;
    
      for(size_t i = 0; i < nVector; i++){
        // Subvector size
        memcpy(&tmpSize, stream + offset, sizeof(size_t));
        offset += sizeof(size_t);
    
        // Alloc Subvector
        dest[i].resize(tmpSize);
    
        // Populate Subvector
        for(size_t j = 0; j < tmpSize; j++){
          memcpy(&tmpVal, stream + offset, sizeof(size_t));
          offset += sizeof(size_t);
    
          dest[i][j] = tmpVal;
        }
      }
    
      // Return size //
      return offset;
    }
    
    size_t ReferenceSpace::
        unserialize(const char* stream,
                    std::vector<std::vector<std::vector<size_t> > >& dest){
      // Vector size //
      size_t nVector;
      size_t offset = 0;
    
      memcpy(&nVector, stream + offset, sizeof(size_t));
      offset += sizeof(size_t);
    
      // Alloc Vector //
      dest.resize(nVector);
    
      // Populate //
      size_t tmpSize;
      size_t tmpSubSize;
      size_t tmpVal;
    
      for(size_t i = 0; i < nVector; i++){
        // Subvector size
        memcpy(&tmpSize, stream + offset, sizeof(size_t));
        offset += sizeof(size_t);
    
        // Alloc Subvector
        dest[i].resize(tmpSize);
    
        // Populate Subvector
        for(size_t j = 0; j < tmpSize; j++){
          // Subsubvector size
          memcpy(&tmpSubSize, stream + offset, sizeof(size_t));
          offset += sizeof(size_t);
    
          // Alloc Subsubvector
          dest[i][j].resize(tmpSubSize);
    
          // Populate Subsubvector
          for(size_t k = 0; k < tmpSubSize; k++){
            memcpy(&tmpVal, stream + offset,sizeof(size_t));
            offset += sizeof(size_t);
    
            dest[i][j][k] = tmpVal; // AT LAST !! (-:
          }
        }
      }
    
      // Return size //
      return offset;
    }