diff --git a/ConSol/ConSol/Mesh/FastHybridMesh.cpp b/ConSol/ConSol/Mesh/FastHybridMesh.cpp
index 0ddcec3..f657e94 100644
--- a/ConSol/ConSol/Mesh/FastHybridMesh.cpp
+++ b/ConSol/ConSol/Mesh/FastHybridMesh.cpp
@@ -1,3282 +1,3290 @@
 //
 // Created by Fabian Moenkeberg on 2019-10-09.
 //
 
 #include "FastHybridMesh.hpp"
 
 #include "Utility/triangle_dunavant_rule.hpp"
 
 #include <cstdlib>
 #include <iostream>
 #include <iomanip>
 #include <cmath>
 #include <cstring>
 #include <algorithm>
 #include <omp.h>
 #include <gmsh.h>
 #include <stdio.h>
 #include "Utility/Quadrature1D/quadrule.hpp"
 
 # define NODE_NUM 4
 
 using namespace std;
 
 void FastHybridMesh::initialize(BcBase *bc, Configuration config) {
     this->bc = bc;
     int internal_num;
     int bcMap_num;
     int edg_num;
     int element_order;
     string gmsh_filename = config.getGridFile();
     int node_num;
     this->nThreads = config.getNthreads();
 
     if (nThreads == 0){
         nThreads = omp_get_max_threads();
     }
 
     std::cout << "  Hybrid Mesh Generation: Running threads: " << nThreads << " out of " << omp_get_max_threads() << "\n";
 
     cout << "\n";
     cout << "  Load mesh data from a file.\n";
 
     gmsh::initialize();
     gmsh::option::setNumber("General.Terminal", 1);
     gmsh::open(gmsh_filename);
     gmsh::model::mesh::generate(3);
     ndims = config.getDim();
 
     std::vector<std::pair<int, int> > entities;
     gmsh::model::getEntities(entities);
 
     std::vector<int> eleTypes;
     gmsh::model::mesh::getElementTypes(eleTypes, ndims);
     int nEleTypes = eleTypes.size();
 
 //    std::vector<std::pair<int, int> > physicalGroups2D;
 //    gmsh::model::getPhysicalGroups(physicalGroups2D, 2);
 
 
     std::string entityName;
     for (int i = 0; i < nEleTypes; i++) {
         int eleType2D = eleTypes[i];
         std::string name;
         int dim, order, numNodes;
         std::vector<double> paramCoord;
         gmsh::model::mesh::getElementProperties(eleType2D, name, dim, order,
                                                 numNodes, paramCoord);
         gmsh::logger::write("2D elements are of type '" + name + "' (type = " +
                             std::to_string(eleType2D) + ") with N = " + std::to_string(numNodes) + " Nodes .");
 //        std::vector<std::pair<int, int> > entities2;
 //        gmsh::model::getEntities(entities2, 2);
 
         gmsh::model::getEntities(entities, 2);
         for (std::size_t i = 0; i < entities.size(); i++) {
             int s = entities[i].second;
             std::vector<std::size_t> elementTags, nodeTags;
 
             gmsh::model::getPhysicalName(2, entities[i].second, entityName);
             gmsh::model::mesh::getElementsByType(eleType2D, elementTags, nodeTags, entities[i].second);
             gmsh::model::getPhysicalName(entities[i].first, entities[i].second, entityName);
             gmsh::model::mesh::getElementsByType(eleType2D, elementTags, nodeTags, s);
             gmsh::logger::write("- " + std::to_string(elementTags.size()) +
                                 " elements in surface " + std::to_string(s) + " (" + entityName + ")");
 
             // get the nodes on the edges of the 2D elements
             std::vector<std::size_t> nodes;
             gmsh::model::mesh::getElementEdgeNodes(eleType2D, nodes, s);
 
             // create a new discrete entity of dimension 1
             int c = gmsh::model::addDiscreteEntity(1);
 
             // and add new 1D elements to it, for all edges
             int eleType1D = gmsh::model::mesh::getElementType("line", order);
             gmsh::model::mesh::addElementsByType(c, eleType1D, {}, nodes);
 
             // this will create two 1D elements for each edge; to create unique elements
             // it would be useful to call getElementEdgeNodes() with the extra `primary'
             // argument set to 'true' (to only get start/end nodes even in the
             // high-order case, i.e. consider topological edges), then sort them and
             // make them unique.
 
             // this could be enriched with additional info: each topological edge could
             // be associated with the tag of its parent element; in the sorting process
             // (eliminating duplicates) a second tag can be associated for internal
             // edges, allowing to keep track of neighbors
         }
     }
     std::vector<std::pair<int, int> > physicalGroups;
     gmsh::model::getPhysicalGroups(physicalGroups, 1);
     std::string phyGroupName;
     std::vector<int> physicalTags;
 
     std::vector<std::size_t> nodeTags;
     std::vector<double> coord;
     std::vector<double> parametricCoord;
     gmsh::model::mesh::getNodes(nodeTags, coord, parametricCoord);
 
     ngc = floor((config.getOrder()-1)/2.0);
 
     ngc = 2;
 
     int nNodes = coord.size()/3;
     xCoord.resize(nNodes,std::vector<double>(2));
 
     for (int i = 0; i < nNodes; i++){
         xCoord[nodeTags[i]-1][0] = coord[3*i];
         xCoord[nodeTags[i]-1][1] = coord[3*i+1];
     }
 
     std::vector<std::size_t> nodeTagsP;
     std::vector<double> coordP;
     std::string physicalName;
     std::vector<bcName> listPhysicalNameEdges;
     std::vector<std::vector<std::size_t>> listPhysicalNodeTagsEdges;
     listPhysicalNameEdges.resize(0);
     listPhysicalNodeTagsEdges.resize(0,std::vector<std::size_t>(0));
     int nPhysicalGroups = physicalGroups.size();
     for (int i = 0; i < nPhysicalGroups; i++){
         gmsh::model::mesh::getNodesForPhysicalGroup(physicalGroups[i].first, physicalGroups[i].second, nodeTagsP,
                                                     coordP);
         gmsh::model::getPhysicalName(physicalGroups[i].first, physicalGroups[i].second, physicalName);
         if (physicalName == "inflow") {
             listPhysicalNameEdges.push_back(INFLOW);
         }else if (physicalName == "outflow"){
             listPhysicalNameEdges.push_back(OUTFLOW);
         }else if (physicalName == "slipWall"){
             listPhysicalNameEdges.push_back(SLIPWALL);
         }else if (physicalName == "exact"){
             listPhysicalNameEdges.push_back(EXACT);
         }else if (physicalName == "periodic"){
             listPhysicalNameEdges.push_back(PERIODIC);
         }else if (physicalName == "special"){
             listPhysicalNameEdges.push_back(SPECIAL);
         }else if (physicalName == "q2i"){
             listPhysicalNameEdges.push_back(Q2I);
         }else if (physicalName == "t2i"){
             listPhysicalNameEdges.push_back(T2I);
         }else if (physicalName == "absorb"){
             listPhysicalNameEdges.push_back(ABSORBINGBC);
         }else if (physicalName == "q2q"){
             listPhysicalNameEdges.push_back(Q2Q);
         }else{
             assert(false && "Read FastHybridMesh: unknown BC ");
         }
 
         for (int j = 0; j < nodeTagsP.size(); j++){
             nodeTagsP[j]--;
         }
         listPhysicalNodeTagsEdges.push_back(nodeTagsP);
         gmsh::logger::write("Physical Group of type '"+ physicalName + " with dimension " + std::to_string(physicalGroups[i].first) + ", " + std::to_string(physicalGroups[i].second) + " has N = " +
                             std::to_string(nodeTagsP.size()) + " tags. ");
     }
 
     gmsh::model::getPhysicalGroups(physicalGroups, 2);
     std::vector<MeshType> listPhysicalNameElements;
     std::vector<int> listPhysicalGroupNumber;
     std::vector<std::vector<std::size_t>> listPhysicalNodeTagsElements;
     listPhysicalNameElements.resize(0);
     listPhysicalNodeTagsElements.resize(0,std::vector<std::size_t>(0));
     nPhysicalGroups = physicalGroups.size();
     for (int i = 0; i < nPhysicalGroups; i++){
         gmsh::model::mesh::getNodesForPhysicalGroup(physicalGroups[i].first, physicalGroups[i].second, nodeTagsP,
                                                     coordP);
         gmsh::model::getPhysicalName(physicalGroups[i].first, physicalGroups[i].second, physicalName);
 
         if (physicalName.compare(0,3,"tri") == 0) {
             listPhysicalNameElements.push_back(TRI);
         }else if (physicalName.compare(0,3,"qdr") == 0){
             listPhysicalNameElements.push_back(QUAD);
         }else if (physicalName.compare(0,3,"mid") == 0){
             listPhysicalNameElements.push_back(TRIQUAD);
         }else if (physicalName.compare(0,3,"q2q") == 0){
             listPhysicalNameElements.push_back(QUADQUAD);
         }else if (physicalName.compare(0,3,"rst") == 0){
             listPhysicalNameElements.push_back(RT);
         }else if (physicalName.compare(0,3,"rsq") == 0){
             listPhysicalNameElements.push_back(RQ);
         }else{
             assert(false && "Read FastHybridMesh: unknown BC");
         }
 
         for (int j = 0; j < nodeTagsP.size(); j++){
             nodeTagsP[j]--;
         }
         listPhysicalNodeTagsElements.push_back(nodeTagsP);
         gmsh::logger::write("Physical Group of type '"+ physicalName + " with dimension " + std::to_string(physicalGroups[i].first) + ", " + std::to_string(physicalGroups[i].second) + " has N = " +
                             std::to_string(nodeTagsP.size()) + " tags. ");
         physicalName.erase(0,4);
         listPhysicalGroupNumber.push_back(std::stoi(physicalName));
     }
 
     std::vector<int> elementTypes;
     std::vector<std::vector<std::size_t> > elementTags;
     std::vector<std::vector<std::size_t> > nodeTagsElements;
     gmsh::model::mesh::getElements(elementTypes,elementTags,nodeTagsElements,2);
 
     int nElementTypes = elementTypes.size();
     std::vector<int> elementNumbers(nElementTypes);
     std::vector<int> elementSizes(nElementTypes);
     std::vector<int> nnodes(0);
     edg.resize(0,std::vector<int>(0));
     e2c.resize(0,std::vector<int>(0));
     c2e.resize(0, std::vector<int>(0));
     cells.resize(0, std::vector<int>(0));
     neighbours.resize(0, std::vector<int>(0));
     nnodes.resize(0);
 
     int nFaces = 0;
     int indFace = 0;
     int nElements0 = 0;
     int iElement = -1;
 
     std::vector<std::size_t> elementTagsByType;
     std::vector<std::size_t > nodeTagsElementsByType;
     submeshCell0.resize(1,0);
 
     std::vector<MeshType> listPhysicalNameElements2(0);
     std::vector<int> listPhysicalGroupNumber2(0);
     std::vector<int> pt;
     int nsubdomains ;
     for (int i = 0; i < nElementTypes; i++) {
         for ( int jj = 0; jj < listPhysicalNodeTagsElements.size(); jj++) {
 //            if (jj == 12){
 //                std::cout <<jj;
 //            }
             int eleType2D = elementTypes[i];
             std::string name;
             int dim, order, numNodes;
             std::vector<double> paramCoord;
             gmsh::model::mesh::getElementProperties(eleType2D, name, dim, order,
                                                     numNodes, paramCoord);
             gmsh::model::getEntitiesForPhysicalGroup(2, jj+1, pt);
             nsubdomains = pt.size();
             nElements0 = 0;
             for (int isD = 0; isD < nsubdomains; isD++) {
                 elementTagsByType.resize(0);
                 nodeTagsElementsByType.resize(0);
                 gmsh::model::mesh::getElementsByType(eleType2D, elementTagsByType, nodeTagsElementsByType, pt[isD]);
 
                 gmsh::model::getPhysicalName(2, jj + 1, entityName);
                 int nElements = elementTagsByType.size();
                 gmsh::logger::write("2D elements are of type '" + name + "' (type = " +
                                     std::to_string(eleType2D) + ") with N = " + std::to_string(numNodes) +
                                     " Nodes, we have " +
                                     std::to_string(nElements) + " elements.");
 
 //            std::vector<int> map(2 * numNodes);
 //            for (int ii = 0; ii < numNodes; ii++) {
 //                map[2 * ii] = ii;
 //                map[2 * ii + 1] = ii + 1;
 //            }
 //            map[2 * numNodes - 1] = 0;
                 std::vector<int> cell0(numNodes);
                 std::vector<int> face(2);
 
                 bool faceExists = 0;
                 elementNumbers[i] = nElements;
                 elementSizes[i] = numNodes;
 
 
                 for (int ii = 0; ii < nElements; ii++) {
                     iElement++;
 
                     nnodes.push_back(numNodes);
 
                     for (int j = 0; j < numNodes; j++) {
                         cell0[j] = nodeTagsElementsByType[numNodes * ii + j] - 1;
                     }
                     cells.push_back(cell0);
                 }
 
                 nElements0 += nElements;
             }
             if (nElements0 > 0){
                 submeshCell0.push_back(nElements0+submeshCell0.back());
                 listPhysicalNameElements2.push_back(listPhysicalNameElements[jj]);
                 listPhysicalGroupNumber2.push_back(listPhysicalGroupNumber[jj]);
             }
         }
     }
 
 // Sort cells such that first all TRI than all QUAD than all TRIQUAD, R0, R1 and finally all QUADQUADs
     gmsh::logger::write("Sort cells: TRI first, QUADs middle, TRIQUADs in the end.");
 
     std::vector<int> iTRI(0);
     std::vector<int> iQUAD(0);
     std::vector<int> iTRIQUAD(0);
     std::vector<int> iQUADQUAD(0);
     std::vector<int> iRQ(0);
     std::vector<int> iRT(0);
     std::vector<std::vector<std::vector<int>>> cellsSplit(0, std::vector<std::vector<int>>(0, std::vector<int>(0)));
     std::vector<int> listPhysicalGroupNumber_copy(listPhysicalGroupNumber);
     for (int ii = 0; ii < listPhysicalNameElements.size(); ii++){
         cellsSplit.push_back(std::vector<std::vector<int>>(cells.begin() + submeshCell0[ii], cells.begin() + submeshCell0[ii+1]));
 
         if (listPhysicalNameElements[ii] == TRI){
             iTRI.push_back(ii);
         }else if (listPhysicalNameElements[ii] == QUAD){
             iQUAD.push_back(ii);
         }else if (listPhysicalNameElements[ii] == TRIQUAD){
             iTRIQUAD.push_back(ii);
         }else if (listPhysicalNameElements[ii] == QUADQUAD){
             iQUADQUAD.push_back(ii);
         }else if (listPhysicalNameElements[ii] == RT){
             iRT.push_back(ii);
         }else if (listPhysicalNameElements[ii] == RQ){
             iRQ.push_back(ii);
         }
     }
 
     nDomainsTRI = iTRI.size();
     nDomainsQUAD = iQUAD.size();
     nDomainsTRIQUAD = iTRIQUAD.size();
     nDomainsQUADQUAD = iQUADQUAD.size();
     nDomainsRQ = iRQ.size();
     nDomainsRT = iRT.size();
 
     cells.resize(0, std::vector<int>(0));
     listPhysicalNameElements.resize(0);
     listPhysicalGroupNumber.resize(0);
     submeshCell0.resize(1,0);
     int nn = 0;
     for (int ii = 0 ; ii < nDomainsTRI; ii++){
         cells.insert(cells.begin()+nn, cellsSplit[iTRI[ii]].begin(), cellsSplit[iTRI[ii]].end());
         nn += cellsSplit[iTRI[ii]].size();
         listPhysicalNameElements.push_back(TRI);
         submeshCell0.push_back(nn);
         listPhysicalGroupNumber.push_back(listPhysicalGroupNumber_copy[iTRI[ii]]);
     }
 
     for (int ii = 0 ; ii < nDomainsQUAD; ii++){
         cells.insert(cells.begin()+nn, cellsSplit[iQUAD[ii]].begin(), cellsSplit[iQUAD[ii]].end());
         nn += cellsSplit[iQUAD[ii]].size();
         listPhysicalNameElements.push_back(QUAD);
         submeshCell0.push_back(nn);
         listPhysicalGroupNumber.push_back(listPhysicalGroupNumber_copy[iQUAD[ii]]);
     }
 
     for (int ii = 0 ; ii < nDomainsTRIQUAD; ii++){
         cells.insert(cells.begin()+nn, cellsSplit[iTRIQUAD[ii]].begin(), cellsSplit[iTRIQUAD[ii]].end());
         nn += cellsSplit[iTRIQUAD[ii]].size();
         listPhysicalNameElements.push_back(TRIQUAD);
         submeshCell0.push_back(nn);
         listPhysicalGroupNumber.push_back(listPhysicalGroupNumber_copy[iTRIQUAD[ii]]);
     }
 
     for (int ii = 0 ; ii < nDomainsRT; ii++){
         cells.insert(cells.begin()+nn, cellsSplit[iRT[ii]].begin(), cellsSplit[iRT[ii]].end());
         nn += cellsSplit[iRT[ii]].size();
         listPhysicalNameElements.push_back(RT);
         submeshCell0.push_back(nn);
         listPhysicalGroupNumber.push_back(listPhysicalGroupNumber_copy[iRT[ii]]);
     }
 
     for (int ii = 0 ; ii < nDomainsQUADQUAD; ii++){
         cells.insert(cells.begin()+nn, cellsSplit[iQUADQUAD[ii]].begin(), cellsSplit[iQUADQUAD[ii]].end());
         nn += cellsSplit[iQUADQUAD[ii]].size();
         listPhysicalNameElements.push_back(QUADQUAD);
         submeshCell0.push_back(nn);
         listPhysicalGroupNumber.push_back(listPhysicalGroupNumber_copy[iQUADQUAD[ii]]);
     }
 
 
     for (int ii = 0 ; ii < nDomainsRQ; ii++){
         cells.insert(cells.begin()+nn, cellsSplit[iRQ[ii]].begin(), cellsSplit[iRQ[ii]].end());
         nn += cellsSplit[iRQ[ii]].size();
         listPhysicalNameElements.push_back(RQ);
         submeshCell0.push_back(nn);
         listPhysicalGroupNumber.push_back(listPhysicalGroupNumber_copy[iRQ[ii]]);
     }
 
 // convert QUAD domains into structured grid.
     gmsh::logger::write("Convert QUAD domains into structured grid.");
     std::vector<int> TRIQUAD2QUAD(nDomainsTRIQUAD);
     std::vector<std::vector<int>> QUADQUAD2QUAD(nDomainsQUADQUAD, std::vector<int>(0));
     std::vector<std::vector<int>> RT2QUADQUAD(nDomainsRT, std::vector<int>(0));
     std::vector<std::vector<int>> RT2QUAD(nDomainsRT, std::vector<int>(0));
     std::vector<std::vector<int>> RQ2QUADQUAD(nDomainsRQ, std::vector<int>(0));
 //    std::vector<std::vector<int>> coordMap(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ, std::vector<int>(0));
     std::vector<int>coordmap(0);
     std::vector<int>listCoordmap(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ + 1,0);
     double eps0 = 1;
     convertQUAD2STRUCTURED(TRIQUAD2QUAD, eps0, coordmap, listCoordmap);
 
     bcMap.resize(0);
     edgesAtBC.resize(0, std::vector<int>(0));
     std::vector<std::vector<int>> bcMap_loc(0, std::vector<int> (0));
 
 // Add ghost cells to each QUAD part, by going through the list of each TRIQUAD domain and add it to the parallel QUAD domain.
     gmsh::logger::write("Add ghost cells to each QUAD part, by going through the list of each TRIQUAD domain and add it to the parallel QUAD domain.");
     std::vector<std::vector<std::vector<int>>> cells_add(listPhysicalNameElements.size()-nDomainsQUAD-nDomainsTRI, std::vector<std::vector<int>>(0,std::vector<int>(0)));
     addGhostCells2QUAD(listPhysicalNameElements, cells_add, eps0, TRIQUAD2QUAD, bcMap_loc);
 
     addGhostCells2QUAD_2(listPhysicalNameElements, cells_add, eps0, QUADQUAD2QUAD, bcMap_loc);
 
 // Add ghost cells to each RT
     addGhostCells2RT(listPhysicalNameElements, cells_add, eps0, RT2QUADQUAD, bcMap_loc, QUADQUAD2QUAD, RT2QUAD);
 
 // Add ghost cells to each RQ
     addGhostCells2RQ(listPhysicalNameElements, cells_add, eps0, RQ2QUADQUAD, bcMap_loc);
 
 // Split cells in TRIQUAD and RT
     splitCellsForRBF(cells_add);
 
     // Add cells_add cells to cells
     gmsh::logger::write("Add cells_add cells to cells.");
     for (int i = nDomainsQUAD + nDomainsTRI; i <  nDomainsQUAD + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT;i++){
         cells.insert(cells.begin()+submeshCell0[i+1],cells_add[i - (nDomainsQUAD + nDomainsTRI)].begin(), cells_add[i - (nDomainsQUAD + nDomainsTRI)].end());
         for (int l = i+1; l < submeshCell0.size(); l++){
             submeshCell0[l] += cells_add[i - (nDomainsQUAD + nDomainsTRI)].size();
         }
     }
 
 
     // Remove QUADs, QUADQUADs and RQs from cells
     gmsh::logger::write("Remove QUADs, QUADQUADs and RQs from cells.");
 
     for (int i = 0; i < nDomainsQUADQUAD + nDomainsRQ; i++){
         cells.erase(cells.begin() + submeshCell0[nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT],cells.begin() + (submeshCell0[nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT]
         + submeshCell0[i + nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT + 1] - submeshCell0[i + nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT]));
     }
 
     for (int i = 0; i < nDomainsQUAD; i++){
         cells.erase(cells.begin() + submeshCell0[nDomainsTRI], cells.begin() + (submeshCell0[nDomainsTRI] + submeshCell0[i + nDomainsTRI + 1] - submeshCell0[i + nDomainsTRI]));
     }
 
     // Change actual size of QUADQUADs
     int i0 = nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT;
     for (int i = 0; i < nDomainsQUADQUAD; i++){
         submeshCell0[i0 + i + 1] = submeshCell0[i0 + i] + (nX[i + nDomainsQUAD] + 2*ngc)*(nY[i + nDomainsQUAD] + 2*ngc);
     }
 
     // Change actual size of RQ
     i0 = nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT + nDomainsQUADQUAD;
     for (int i = 0; i < nDomainsRQ; i++){
         submeshCell0[i0 + i + 1] = submeshCell0[i0 + i] + (nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 2*ngc)*(nY[i + nDomainsQUAD + nDomainsQUADQUAD] + 2*ngc);
     }
 
     // change order of domains such that QUADS are at the end
     gmsh::logger::write("Change order of domains such that QUADS are at the end.");
     std::vector<int> diff(0);
     for (int i = nDomainsTRI; i <  submeshCell0.size()-1; i++){
         diff.push_back((submeshCell0[i+1]-submeshCell0[i]));
     }
     for (int i = 0; i < nDomainsTRIQUAD + nDomainsRT + nDomainsRQ + nDomainsQUADQUAD; i++){
         submeshCell0[i + nDomainsTRI+1] = submeshCell0[i + nDomainsTRI] + diff[nDomainsQUAD + i];
     }
     for (int i = 0; i < nDomainsQUAD; i++){
         submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsRQ + nDomainsQUADQUAD + 1] = submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD
                      + nDomainsRT + nDomainsRQ + nDomainsQUADQUAD] + (nX[i]+2*ngc)*(nY[i]+2*ngc);
     }
 
     for (int i = 0; i < nDomainsQUAD; i++){
         nX.push_back(nX[i]);
         nY.push_back(nY[i]);
         dX.push_back(dX[i]);
         dY.push_back(dY[i]);
         xMin.push_back(xMin[i]);
         xMax.push_back(xMax[i]);
         yMin.push_back(yMin[i]);
         yMax.push_back(yMax[i]);
 //        coordMap.push_back(coordMap[i]);
     }
 
     nX.erase(nX.begin(), nX.begin()+nDomainsQUAD);
     nY.erase(nY.begin(), nY.begin()+nDomainsQUAD);
     dX.erase(dX.begin(), dX.begin()+nDomainsQUAD);
     dY.erase(dY.begin(), dY.begin()+nDomainsQUAD);
     xMin.erase(xMin.begin(), xMin.begin()+nDomainsQUAD);
     xMax.erase(xMax.begin(), xMax.begin()+nDomainsQUAD);
     yMin.erase(yMin.begin(), yMin.begin()+nDomainsQUAD);
     yMax.erase(yMax.begin(), yMax.begin()+nDomainsQUAD);
 //    coordMap.erase(coordMap.begin(), coordMap.begin()+nDomainsQUAD);
 
 
     // set global bcMap
     gmsh::logger::write("Set global bcMap.");
     setGlobalBcMap(bcMap_loc, TRIQUAD2QUAD, QUADQUAD2QUAD, RQ2QUADQUAD, RT2QUADQUAD, RT2QUAD);
 
     // find double cells
     int nc = cells.size();
     std::vector<int>cellMap(2*nc,-1);
     for (int i = 0; i < nc; i++){
         for (int j = i+1; j < nc; j++){
             if (cells[i][0] == cells[j][0]&&cells[i][1] == cells[j][1]&&cells[i][2] == cells[j][2]){
                 cellMap[2*j] = i;
                 cellMap[2*j + 1] += 1;
             }
         }
     }
 
     for (int i = 0; i < nc; i++) {
         if (cellMap[2*i]!=-1){
             cells[i][0] = 0;
             cells[i][1] = 0;
             cells[i][2] = 0;
         }
     }
 
     int nbm = bcMap.size();
     for (int i = 0; i < nbm; i++){
         if (bcMap[i]!=-1&&bcMap[i]<nc){
             if (cellMap[2*bcMap[i]]!=-1){
                 bcMap[i] = cellMap[2*bcMap[i]];
             }
         }
     }
 
     nbm = bcMap2.size();
     for (int i = 0; i < nbm; i++){
         if (bcMap2[i]!=-1&&bcMap2[i]<nc){
             if (cellMap[2*bcMap2[i]]!=-1){
                 bcMap2[i] = cellMap[2*bcMap2[i]];
             }
         }
     }
 
     // set internal map for TRI and TRIQUAD and RT
     int nCells0 = 0;
     for (int i = 0; i < nDomainsTRI; i++){
         nCells0 += submeshCell0[i+1] - submeshCell0[i];
     }
 
     for (int i = 0; i < nDomainsTRIQUAD; i++){
         int n = submeshCell0[i+1+nDomainsTRI] - submeshCell0[i+nDomainsTRI];
         nCells0 += n/2;
     }
 
     for (int i = 0; i < nDomainsRT; i++){
 //        int n = (submeshCell0[i + 1 + nDomainsTRI + nDomainsTRIQUAD] - submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD])/(RT2QUADQUAD[i].size() + RT2QUAD[i].size() + 1);
         int n = ngc*ngc*2;
         nCells0 += n;
     }
 
 
     nCells = nCells0;
 
     // set internal map for RQ, QUADQUAD, QUAD
     gmsh::logger::write("Set internal map for RQ, QUADQUAD, QUAD.");
     for (int i = 0; i < nDomainsQUADQUAD + nDomainsRQ + nDomainsQUAD; i++){
         nCells0 += nX[i]*nY[i];
     }
 
     internal.resize(nCells0);
     int icell = 0;
     gmsh::logger::write("Set internal map for TRI and TRIQUAD and RT.");
     for (int i = 0; i < nDomainsTRI; i++){
         for (int j = submeshCell0[i]; j < submeshCell0[i+1]; j++){
             internal[icell] = j;
             icell++;
         }
     }
 
     for (int i = 0; i < nDomainsTRIQUAD; i++){
         int n = submeshCell0[i+1+nDomainsTRI] - submeshCell0[i+nDomainsTRI];
         for (int j = 0; j < n/2; j++){
             internal[icell] = j + submeshCell0[i+nDomainsTRI];
             icell++;
         }
     }
 
     for (int i = 0; i < nDomainsRT; i++){
 //        int n = (submeshCell0[i + 1 + nDomainsTRI + nDomainsTRIQUAD] - submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD])/(RT2QUADQUAD[i].size() + RT2QUAD[i].size() + 1);
         int n = ngc*ngc*2;
         for (int j = 0; j < n; j++){
             internal[icell] = j + submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD];
             icell++;
         }
     }
 
 
     // set internal map for RQ, QUADQUAD, QUAD
     gmsh::logger::write("Set internal map for RQ, QUADQUAD, QUAD.");
     for (int i = 0; i < nDomainsQUADQUAD; i++){
         for (int i0 = 0; i0 < nX[i]; i0++){
             for (int j0 = 0; j0 < nY[i]; j0++){
                 internal[icell] = submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT] + i0 + ngc + (j0 + ngc)*(nX[i] + 2*ngc);
                 icell++;
             }
         }
     }
 
     for (int i = 0; i < nDomainsRQ; i++){
         for (int i0 = 0; i0 < nX[i + nDomainsQUADQUAD]; i0++){
             for (int j0 = 0; j0 < nY[i + nDomainsQUADQUAD]; j0++){
                 internal[icell] = submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsQUADQUAD] + i0 + ngc + (j0 + ngc)*(nX[i + nDomainsQUADQUAD] + 2*ngc);
                 icell++;
             }
         }
     }
 
     for (int i = 0; i < nDomainsQUAD; i++){
         for (int i0 = 0; i0 < nX[i + nDomainsRQ + nDomainsQUADQUAD]; i0++){
             for (int j0 = 0; j0 < nY[i + nDomainsRQ + nDomainsQUADQUAD]; j0++){
                 internal[icell] = submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsRQ + nDomainsQUADQUAD] + i0 + ngc + (j0 + ngc)*(nX[i + nDomainsRQ + nDomainsQUADQUAD] + 2*ngc);
                 icell++;
             }
         }
     }
 
 //    int nCells = 0;
 //
 //    gmsh::logger::write("Set internal map for TRI and TRIQUAD and RT.");
 //    for (int i = 0; i < nDomainsTRI; i++){
 //        for (int j = submeshCell0[i]; j < submeshCell0[i+1]; j++){
 //            internal.push_back(j);
 //        }
 //    }
 //
 //    for (int i = 0; i < nDomainsTRIQUAD; i++){
 //        int n = submeshCell0[i+1+nDomainsTRI] - submeshCell0[i+nDomainsTRI];
 //        for (int j = 0; j < n/2; j++){
 //            internal.push_back(j + submeshCell0[i+nDomainsTRI]);
 //        }
 //    }
 //
 //    for (int i = 0; i < nDomainsRT; i++){
 ////        int n = (submeshCell0[i + 1 + nDomainsTRI + nDomainsTRIQUAD] - submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD])/(RT2QUADQUAD[i].size() + RT2QUAD[i].size() + 1);
 //        int n = ngc*ngc*2;
 //        for (int j = 0; j < n; j++){
 //            internal.push_back(j + submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD]);
 //        }
 //    }
 //
 //    nCells = internal.size();
 //
 //
 //    // set internal map for RQ, QUADQUAD, QUAD
 //    gmsh::logger::write("Set internal map for RQ, QUADQUAD, QUAD.");
 //    for (int i = 0; i < nDomainsQUADQUAD; i++){
 //        for (int i0 = 0; i0 < nX[i]; i0++){
 //            for (int j0 = 0; j0 < nY[i]; j0++){
 //                internal.push_back(submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT] + i0 + ngc + (j0 + ngc)*(nX[i] + 2*ngc));
 //            }
 //        }
 //    }
 //
 //    for (int i = 0; i < nDomainsRQ; i++){
 //        for (int i0 = 0; i0 < nX[i + nDomainsQUADQUAD]; i0++){
 //            for (int j0 = 0; j0 < nY[i + nDomainsQUADQUAD]; j0++){
 //                internal.push_back(submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsQUADQUAD] + i0 + ngc + (j0 + ngc)*(nX[i + nDomainsQUADQUAD] + 2*ngc));
 //            }
 //        }
 //    }
 //
 //    for (int i = 0; i < nDomainsQUAD; i++){
 //        for (int i0 = 0; i0 < nX[i + nDomainsRQ + nDomainsQUADQUAD]; i0++){
 //            for (int j0 = 0; j0 < nY[i + nDomainsRQ + nDomainsQUADQUAD]; j0++){
 //                internal.push_back(submeshCell0[i + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsRQ + nDomainsQUADQUAD] + i0 + ngc + (j0 + ngc)*(nX[i + nDomainsRQ + nDomainsQUADQUAD] + 2*ngc));
 //            }
 //        }
 //    }
 
 //    std::vector<int> notInternal(0);
 //    int iCell = 0;
 //    for (int i = 0; i < nCells; i++){
 //        if(internal[i] > iCell){
 //            for (; iCell < internal[i]; iCell++){
 //                notInternal.push_back(iCell);
 //            }
 //            iCell++;
 //        }else{
 //            iCell++;
 //        }
 //    }
 //
 //    int nElements = cells.size();
 //    for (;iCell < nElements; iCell++){
 //        notInternal.push_back(iCell);
 //    }
 //    iElement = -1;
 //    for (int ii = 0; ii < nCells; ii++) {
 ////        iElement++;
 //        iElement = internal[ii];
 //        int numNodes = cells[ii].size();
 //        std::vector<int> cell0(numNodes);
 //        std::vector<int> face(2);
 //        nnodes.push_back(numNodes);
 //        c2e.push_back(cell0);
 ////        std::vector<int> map(2 * numNodes);
 ////        for (int ii = 0; ii < numNodes; ii++) {
 ////            map[2 * ii] = ii;
 ////            map[2 * ii + 1] = ii + 1;
 ////        }
 ////        map[2 * numNodes - 1] = 0;
 //        bool faceExists = 0;
 //
 //        for (int j = 0; j < numNodes; j++) {
 //            c2e[ii][j] = -1;
 ////            cell0[j] = nodeTagsElementsByType[numNodes * ii + j] - 1;
 ////            face[0] = nodeTagsElementsByType[numNodes * ii + map[2 * j]] - 1;
 ////            face[1] = nodeTagsElementsByType[numNodes * ii + map[2 * j + 1]] - 1;
 //            face[0] = cells[iElement][j];
 //            face[1] = cells[iElement][(j+1)%numNodes];
 //            for (int jj = 0; jj < nFaces; jj++) {
 //                if (((edg[jj][0] == face[0]) && (edg[jj][1] == face[1])) ||
 //                    ((edg[jj][0] == face[1]) && (edg[jj][1] == face[0]))) {
 //                    faceExists = true;
 //                    indFace = jj;
 //                }
 //            }
 //
 //            if (faceExists) {
 //                c2e[ii][j] = indFace;
 //                e2c[indFace][1] = iElement;
 //            } else {
 //                nFaces++;
 //                c2e[ii][j] = nFaces - 1;
 //                edg.push_back(face);
 //                e2c.push_back({iElement, -1});
 //            }
 //            faceExists = 0;
 //        }
 //    }
 //    int nCellsNot = notInternal.size();
 //    for (int ii = 0; ii < nCellsNot; ii++) {
 ////        iElement++;
 //        iElement = notInternal[ii];
 //        int numNodes = cells[ii].size();
 //        std::vector<int> cell0(numNodes);
 //        std::vector<int> face(2);
 //        nnodes.push_back(numNodes);
 //        c2e.push_back(cell0);
 ////        std::vector<int> map(2 * numNodes);
 ////        for (int ii = 0; ii < numNodes; ii++) {
 ////            map[2 * ii] = ii;
 ////            map[2 * ii + 1] = ii + 1;
 ////        }
 ////        map[2 * numNodes - 1] = 0;
 //        bool faceExists = 0;
 //
 //        for (int j = 0; j < numNodes; j++) {
 //            c2e[ii + nCells][j] = -1;
 ////            cell0[j] = nodeTagsElementsByType[numNodes * ii + j] - 1;
 ////            face[0] = nodeTagsElementsByType[numNodes * ii + map[2 * j]] - 1;
 ////            face[1] = nodeTagsElementsByType[numNodes * ii + map[2 * j + 1]] - 1;
 //            face[0] = cells[iElement][j];
 //            face[1] = cells[iElement][(j+1)%numNodes];
 //            for (int jj = 0; jj < nFaces; jj++) {
 //                if (((edg[jj][0] == face[0]) && (edg[jj][1] == face[1])) ||
 //                    ((edg[jj][0] == face[1]) && (edg[jj][1] == face[0]))) {
 //                    faceExists = true;
 //                    indFace = jj;
 //                }
 //            }
 //
 //            if (faceExists) {
 //                c2e[ii + nCells][j] = indFace;
 //                e2c[indFace][1] = iElement;
 //            } else {
 //                nFaces++;
 //                c2e[ii + nCells][j] = nFaces - 1;
 //                edg.push_back(face);
 //                e2c.push_back({iElement, -1});
 //            }
 //            faceExists = 0;
 //        }
 //    }
 
 
     // define edg, c2e, e2c mapping for TRI/TRIQUAD
 
     int nElements = cells.size();
     std::vector<int> isInternal(nElements);
     for (int i = 0; i < nCells; i++){
         isInternal[internal[i]] = 1;
     }
     iElement = -1;
     for (int ii = 0; ii < nElements; ii++) {
         iElement++;
 //        if (iElement == 105){
 //            std::cout << iElement ;
 //        }
         int numNodes = cells[ii].size();
         std::vector<int> cell0(numNodes);
         std::vector<int> face(2);
         nnodes.push_back(numNodes);
         c2e.push_back(cell0);
 //        std::vector<int> map(2 * numNodes);
 //        for (int ii = 0; ii < numNodes; ii++) {
 //            map[2 * ii] = ii;
 //            map[2 * ii + 1] = ii + 1;
 //        }
 //        map[2 * numNodes - 1] = 0;
         bool faceExists = 0;
 
         for (int j = 0; j < numNodes; j++) {
             c2e[iElement][j] = -1;
 //            cell0[j] = nodeTagsElementsByType[numNodes * ii + j] - 1;
 //            face[0] = nodeTagsElementsByType[numNodes * ii + map[2 * j]] - 1;
 //            face[1] = nodeTagsElementsByType[numNodes * ii + map[2 * j + 1]] - 1;
             face[0] = cells[iElement][j];
             face[1] = cells[iElement][(j+1)%numNodes];
             for (int jj = 0; jj < nFaces; jj++) {
                 if (((edg[jj][0] == face[0]) && (edg[jj][1] == face[1])) ||
                     ((edg[jj][0] == face[1]) && (edg[jj][1] == face[0]))) {
                     faceExists = true;
                     indFace = jj;
                 }
             }
 
             if (faceExists) {
                 c2e[iElement][j] = indFace;
                 if (isInternal[iElement] == 0){
                     e2c[indFace][1] = iElement;
                 }else{
                     e2c[indFace][1] = e2c[indFace][0];
                     e2c[indFace][0] = iElement;
                 }
             } else {
                 nFaces++;
                 c2e[iElement][j] = nFaces - 1;
                 edg.push_back(face);
                 e2c.push_back({iElement, -1});
             }
             faceExists = 0;
         }
     }
 
 
 
 
     // Make sure that the cells are oriented counter-clockwise.
     gmsh::logger::write("Make sure that the cells are oriented counter-clockwise.");
     int numNodes;
     double cp;
     int nElementsTot = cells.size();
     for (int i = 0; i < nElementsTot; i++){
         numNodes = cells[i].size();
         for (int j = 0; j < numNodes; j++){
             cp = (xCoord[cells[i][(j+1)%numNodes]][0] - xCoord[cells[i][j]][0])*(xCoord[cells[i][(j+2)%numNodes]][1] - xCoord[cells[i][(j+1)%numNodes]][1]) -
                  (xCoord[cells[i][(j+1)%numNodes]][1] - xCoord[cells[i][j]][1])*(xCoord[cells[i][(j+2)%numNodes]][0] - xCoord[cells[i][(j+1)%numNodes]][0]);
             if (cp < 0){
                 assert(false && "FastHybridMesh: Wrong orientation of cell in grid !!!!");
             }
         }
     }
 
     // Make sure that the cell to the "left" of each edge sees that edge
     // as being oriented counter-clockwise.
     nFaces = e2c.size();
     std::vector<int> left(nFaces);
     std::vector<int> right(nFaces);
     std::vector<int> cell;
 
     int tmp;
     for (int i = 0; i < nFaces; i++){
         left[i] = e2c[i][0];
         numNodes = cells[left[i]].size();
         bool badOriented = true;
 
         for (int j = 0; j < numNodes-1; j++){
             if ((cells[left[i]][j] == edg[i][0])&&(cells[left[i]][j+1] == edg[i][1])){badOriented = false;}
         }
         if ((cells[left[i]][numNodes-1]==edg[i][0])&&(cells[left[i]][0]==edg[i][1])){badOriented = false;}
 
         if (badOriented){
             tmp = edg[i][1];
             edg[i][1] = edg[i][0];
             edg[i][0] = tmp;
         }
     }
 
     // Calculate neighbours mapping or neighbour mapping
     gmsh::logger::write("Calculate neighbours mapping or neighbour mapping.");
     nElements0 = 0;
     int nCellsTRIQUAD = submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT];
     int n;
     iElement = -1;
 
     for (int ii = 0; ii < nCellsTRIQUAD; ii++) {
         iElement++;
         int numNodes = cells[iElement].size();
         std::vector<int> cell0(numNodes);
         neighbours.push_back(cell0);
 
         for (int j = 0; j < numNodes; j++) {
             n = e2c[c2e[iElement][j]][0];
             if (n == iElement) {
                 neighbours[iElement][j] = e2c[c2e[iElement][j]][1];
             } else {
                 neighbours[iElement][j] = n;
             }
             if (neighbours[iElement][j] == -1) {
                 neighbours[iElement][j] = iElement;
             }
         }
     }
 
 
 
     submeshEdg0.resize(0);
     submeshEdg0.push_back(0);
     submeshEdg0.push_back(edg.size());
     for (int i = 0; i < nDomainsQUADQUAD; i++){
         submeshEdg0.push_back(submeshEdg0[i + 1] + (nY[i]+1)*nX[i] + (nX[i]+1)*nY[i]);
     }
     for (int i = 0; i < nDomainsRQ; i++){
         submeshEdg0.push_back(submeshEdg0[i + 1 + nDomainsQUADQUAD] + (nY[i + nDomainsQUADQUAD]+1)*nX[i + nDomainsQUADQUAD] + (nX[i + nDomainsQUADQUAD]+1)*nY[i + nDomainsQUADQUAD]);
     }
     for (int i = 0; i < nDomainsQUAD; i++){
         submeshEdg0.push_back(submeshEdg0[i + 1 + nDomainsRQ + nDomainsQUADQUAD] + (nY[i + nDomainsRQ + nDomainsQUADQUAD]+1)*nX[i + nDomainsRQ + nDomainsQUADQUAD]
         + (nX[i + nDomainsRQ + nDomainsQUADQUAD]+1)*nY[i + nDomainsRQ + nDomainsQUADQUAD]);
     }
 
 
     // Add for each QUAD, QUADQUAD and RQ on boundary its bcMap and edgesAtBC
     gmsh::logger::write("Add for each QUAD on boundary its bcMap and edgesAtBC.");
 //    edgesAtHybridBC.resize(nDomainsQUAD, std::vector<std::vector<int>>(0,std::vector<int>(0)));
     std::vector<int> edgeInfo(4);
     int nBC = listPhysicalNodeTagsEdges.size();
 
     for (int iii = 0; iii < nDomainsRQ + nDomainsQUADQUAD + nDomainsQUAD; iii++){
         // edges (x,y_0)
         for (int i = 0; i < nX[iii]; i++){
             int j = 0;
 
             edgeInfo[0] = i + (j)*(nX[iii]) + submeshEdg0[iii + 1]; // general map
 //            if (edgeInfo[0] == 346){
 //                std::cout << i;
 //            }
 //        edgeInfo[0] = 0; // local map
             for (int jj = 0; jj < nBC; jj++){
                 int n = listPhysicalNodeTagsEdges[jj].size();
                 int count = 0;
                 for (int ii = 0; ii < n; ii++){
                     if ((coordmap[listCoordmap[iii] + i + (j)*(nX[iii] + 1)] == listPhysicalNodeTagsEdges[jj][ii])||
                         (coordmap[listCoordmap[iii] + i + 1 + (j)*(nX[iii] + 1)] == listPhysicalNodeTagsEdges[jj][ii])){
                         count++;
                     }
 //                    if ((coordMap[iii][i + (j)*(nX[iii] + 1)] == listPhysicalNodeTagsEdges[jj][ii])||
 //                    (coordMap[iii][i + 1 + (j)*(nX[iii] + 1)] == listPhysicalNodeTagsEdges[jj][ii])){
 //                        count++;
 //                    }
                 }
                 if (count == 2){
                     edgeInfo[1] = int(listPhysicalNameEdges[jj]);
                 }
             }
 //            edgeInfo[2] = 0;    // means right side of edge is out & normal is OUTWARDS pointing.
             edgeInfo[2] = 1;    // 1 would mean left side of edge is out & normal is INWARDS pointing.
 //            if (iii < nDomainsRQ + nDomainsQUADQUAD){ //for QUADQUADs and RQs.
 //                edgeInfo[2] = 0;    // means right side of edge is out & normal is OUTWARDS pointing.
 //            }
 //            edgeInfo[3] = i + nY[iii]*(nX[iii]); // general map
 //        edgeInfo[3] = nY; // local map
             edgeInfo[3] = -1; //means for periodic bc not working until set properly...
             if (edgeInfo[1]!=Q2I && edgeInfo[1]!= T2I && edgeInfo[1]!= Q2Q){
                 edgesAtBC.push_back(edgeInfo);
                 for (int jcell = 0; jcell < ngc; jcell++){
                     std::vector<int> bm(3);
 //                    bm[0] = i + ngc + (jcell)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + iii];
 //                    bm[1] = i + ngc + (2*ngc - 1 - jcell)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + iii];
 //                    bm[2] = edgesAtBC.size()-1;
                     bcMap.push_back(i + ngc + (jcell)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     if (edgeInfo[1] == ABSORBINGBC){
                         bcMap.push_back(i + ngc + (ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     }else{
                         bcMap.push_back(i + ngc + (2*ngc - 1 - jcell)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     }
 
                     bcMap.push_back(edgesAtBC.size()-1);
                 }
             }else if (edgeInfo[1]==Q2Q){
                 for (int iii2 = iii+1; iii2 < nDomainsRQ + nDomainsQUADQUAD + nDomainsQUAD; iii2++){
                     if (iii2!=iii){
                         if (std::fabs(xMin[iii] - xMin[iii2])<eps0&&std::fabs(xMax[iii] - xMax[iii2])<eps0&&std::fabs(yMin[iii] + 2*ngc*dY[iii] - yMax[iii2])<eps0){
                             edgeInfo[3] = i + (nY[iii2])*(nX[iii2]) + submeshEdg0[iii2 + 1];
                             edgesAtBC.push_back(edgeInfo);
                             edgeInfo[3] = edgeInfo[0];
                             edgeInfo[0] = i + (nY[iii2])*(nX[iii2]) + submeshEdg0[iii2 + 1];
 //                            edgeInfo[2] = (edgeInfo[2] + 1)%2;
                             edgesAtBC.push_back(edgeInfo);
                             break;
                         }
                     }
                 }
             }else{
                 int nEdg_temp = edg.size();
                 for (int iEdg = 0; iEdg < nEdg_temp; iEdg++){
                     if(std::fabs((i + 0.5 + ngc)*dX[iii] + xMin[iii] - (xCoord[edg[iEdg][0]][0]+xCoord[edg[iEdg][1]][0])/2)<eps0 &&
                     std::fabs((j + ngc)*dY[iii] + yMin[iii] - (xCoord[edg[iEdg][0]][1]+xCoord[edg[iEdg][1]][1])/2)<eps0){
                         edgeInfo[3] = iEdg;
                         edgesAtBC.push_back(edgeInfo);
                         edgeInfo[3] = edgeInfo[0];
                         edgeInfo[0] = iEdg;
                         edgeInfo[2] = (edgeInfo[2] + 1)%2;
                         edgesAtBC.push_back(edgeInfo);
                         break;
                     }
                 }
             }
         }
 
         // edges (x,y_end)
         for (int i = 0; i < nX[iii]; i++){
             int j = nY[iii];
 
             edgeInfo[0] = i + (j)*(nX[iii]) + submeshEdg0[iii + 1]; // general map
 //            if (edgeInfo[0] == 2007){
 //                std::cout << i;
 //            }
 //        edgeInfo[0] = j; // local map
             for (int jj = 0; jj < nBC; jj++){
                 int n = listPhysicalNodeTagsEdges[jj].size();
                 int count = 0;
                 for (int ii = 0; ii < n; ii++){
                     if ((coordmap[listCoordmap[iii] + i+j*(nX[iii]+1)] == listPhysicalNodeTagsEdges[jj][ii])||
                     (coordmap[listCoordmap[iii] + i+1+j*(nX[iii]+1)] == listPhysicalNodeTagsEdges[jj][ii])){
                         count++;
                     }
                 }
                 if (count == 2){
                     edgeInfo[1] = int(listPhysicalNameEdges[jj]);}
             }
             edgeInfo[2] = 1;    // means left side of edge is out & normal is INWARDS pointing.
 //        edgeInfo[2] = 0;    // means right side of edge is out & normal is OUTWARDS pointing.
 //            if (iii < nDomainsRQ + nDomainsQUADQUAD){ //for QUADQUADs and RQs.
 //                edgeInfo[2] = 0;    // means right side of edge is out & normal is OUTWARDS pointing.
 //            }
 //            edgeInfo[3] = i + 0*(nX); //general map
 //        edgeInfo[3] = 0; //local map
             edgeInfo[3] = -1; //means for periodic bc not working until set properly...
             if (edgeInfo[1]!=Q2I && edgeInfo[1]!= T2I && edgeInfo[1]!= Q2Q){
                 edgesAtBC.push_back(edgeInfo);
                 for (int jcell = 0; jcell < ngc; jcell++){
                     std::vector<int> bm(3);
                     bcMap.push_back(i + ngc + (jcell + nY[iii] + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     if (edgeInfo[1] == ABSORBINGBC){
                         bcMap.push_back(i + ngc + ( - 1 + nY[iii] + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     } else{
                         bcMap.push_back(i + ngc + (-jcell - 1 + nY[iii] + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     }
 
                     bcMap.push_back(edgesAtBC.size()-1);
                 }
             }else if (edgeInfo[1]==Q2Q){
                 for (int iii2 = iii+1; iii2 < nDomainsRQ + nDomainsQUADQUAD + nDomainsQUAD; iii2++){
                     if (iii2!=iii){
                         if (std::fabs(xMin[iii] - xMin[iii2])<eps0&&std::fabs(xMax[iii] - xMax[iii2])<eps0&&std::fabs(yMin[iii2] + 2*ngc*dY[iii] - yMax[iii])<eps0){
                             edgeInfo[3] = i + (0)*(nX[iii2]) + submeshEdg0[iii2 + 1];
                             edgesAtBC.push_back(edgeInfo);
                             edgeInfo[3] = edgeInfo[0];
                             edgeInfo[0] = i + (0)*(nX[iii2]) + submeshEdg0[iii2 + 1];
 //                            edgeInfo[2] = (edgeInfo[2] + 1)%2;
                             edgesAtBC.push_back(edgeInfo);
                             break;
                         }
                     }
                 }
             }else{
                 int nEdg_temp = edg.size();
                 for (int iEdg = 0; iEdg < nEdg_temp; iEdg++){
                     if(std::fabs((i + 0.5 + ngc)*dX[iii] + xMin[iii] - (xCoord[edg[iEdg][0]][0]+xCoord[edg[iEdg][1]][0])/2)<eps0 &&
                        std::fabs((j + ngc)*dY[iii] + yMin[iii] - (xCoord[edg[iEdg][0]][1]+xCoord[edg[iEdg][1]][1])/2)<eps0){
                         edgeInfo[3] = iEdg;
                         edgesAtBC.push_back(edgeInfo);
                         edgeInfo[3] = edgeInfo[0];
                         edgeInfo[0] = iEdg;
                         edgeInfo[2] = (edgeInfo[2] + 1)%2;
                         edgesAtBC.push_back(edgeInfo);
                         break;
                     }
                 }
             }
         }
 
         // edges (x_0,y)
         for (int j = 0; j < nY[iii]; j++){
             int i = 0;
 
             edgeInfo[0] = i + (j)*(nX[iii] + 1) + nX[iii]*(nY[iii]+1) + submeshEdg0[iii + 1]; // general map
 //            if (edgeInfo[0] == 2007){
 //                std::cout << i;
 //            }
 //        edgeInfo[0] = i; // local map
             for (int jj = 0; jj < nBC; jj++){
                 int n = listPhysicalNodeTagsEdges[jj].size();
                 int count = 0;
                 for (int ii = 0; ii < n; ii++){
                     if ((coordmap[listCoordmap[iii] + i+j*(nX[iii]+1)] == listPhysicalNodeTagsEdges[jj][ii])|| (coordmap[listCoordmap[iii] + i +(j+1)*(nX[iii]+1)] == listPhysicalNodeTagsEdges[jj][ii])){
                         count++;
                     }
                 }
                 if (count == 2){
                     edgeInfo[1] = int(listPhysicalNameEdges[jj]);}
             }
 //            edgeInfo[2] = 0;    // means right side of edge is out & normal is OUTWARDS pointing.
             edgeInfo[2] = 1;    // 1 would mean left side of edge is out & normal is INWARDS pointing.
 //            if (iii < nDomainsRQ + nDomainsQUADQUAD){ //for QUADQUADs and RQs.
 //                edgeInfo[2] = 0;    // means right side of edge is out & normal is OUTWARDS pointing.
 //            }
 //            edgeInfo[3] = nX + j*(nX+1)+ nX*(nY+1); // general map
 //        edgeInfo[3] = nX; // local map
             edgeInfo[3] = -1; //means for periodic bc not working until set properly...
             if (edgeInfo[1]!=Q2I && edgeInfo[1]!= T2I && edgeInfo[1]!= Q2Q){
                 edgesAtBC.push_back(edgeInfo);
                 for (int icell = 0; icell < ngc; icell++){
                     std::vector<int> bm(3);
                     bcMap.push_back(icell + (j + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     if (edgeInfo[1] == ABSORBINGBC){
                         bcMap.push_back(ngc + (j + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     }else{
                         bcMap.push_back(2*ngc - icell -1 + (j + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     }
                     bcMap.push_back(edgesAtBC.size()-1);
                 }
             }else if (edgeInfo[1]==Q2Q){
                 for (int iii2 = iii+1; iii2 < nDomainsRQ + nDomainsQUADQUAD + nDomainsQUAD; iii2++){
                     if (iii2!=iii){
                         if (std::fabs(yMin[iii] - yMin[iii2])<eps0&&std::fabs(yMax[iii] - yMax[iii2])<eps0&&std::fabs(xMin[iii] + 2*ngc*dX[iii] - xMax[iii2])<eps0){
                             edgeInfo[3] = nX[iii2] + (j)*(nX[iii2] + 1) + nX[iii2]*(nY[iii2]+1) + submeshEdg0[iii2 + 1];
                             edgesAtBC.push_back(edgeInfo);
                             edgeInfo[3] = edgeInfo[0];
                             edgeInfo[0] = nX[iii2] + (j)*(nX[iii2] + 1) + nX[iii2]*(nY[iii2]+1) + submeshEdg0[iii2 + 1];
 //                            edgeInfo[2] = (edgeInfo[2] + 1)%2;
                             edgesAtBC.push_back(edgeInfo);
                             break;
                         }
                     }
                 }
             }else{
                 int nEdg_temp = edg.size();
                 for (int iEdg = 0; iEdg < nEdg_temp; iEdg++){
                     if(std::fabs((i + ngc)*dX[iii] + xMin[iii] - (xCoord[edg[iEdg][0]][0]+xCoord[edg[iEdg][1]][0])/2)<eps0 &&
                        std::fabs((j + 0.5 + ngc)*dY[iii] + yMin[iii] - (xCoord[edg[iEdg][0]][1]+xCoord[edg[iEdg][1]][1])/2)<eps0){
                         edgeInfo[3] = iEdg;
                         edgesAtBC.push_back(edgeInfo);
                         edgeInfo[3] = edgeInfo[0];
                         edgeInfo[0] = iEdg;
                         edgeInfo[2] = (edgeInfo[2] + 1)%2;
                         edgesAtBC.push_back(edgeInfo);
                         break;
                     }
                 }
             }
         }
 
         // edges (x_end,y)
         for (int j = 0; j < nY[iii]; j++){
             int i = nX[iii];
 
             edgeInfo[0] = i + (j)*(nX[iii] + 1) + nX[iii]*(nY[iii]+1) + submeshEdg0[iii + 1]; // general map
 //            if (edgeInfo[0] == 2007){
 //                std::cout << i;
 //            }
 //        edgeInfo[0] = i; // local map
             for (int jj = 0; jj < nBC; jj++){
                 int n = listPhysicalNodeTagsEdges[jj].size();
                 int count = 0;
                 for (int ii = 0; ii < n; ii++){
                     if ((coordmap[listCoordmap[iii] + i+j*(nX[iii]+1)] == listPhysicalNodeTagsEdges[jj][ii])|| (coordmap[listCoordmap[iii] + i+(j+1)*(nX[iii]+1)] == listPhysicalNodeTagsEdges[jj][ii])){
                         count++;
                     }
                 }
                 if (count == 2){
                     edgeInfo[1] = int(listPhysicalNameEdges[jj]);}
             }
             edgeInfo[2] = 1;    // means left side of edge is out & normal is INWARDS pointing.
 //        edgeInfo[2] = 0;    // means right side of edge is out & normal is OUTWARDS pointing.
 //            if (iii < nDomainsRQ + nDomainsQUADQUAD){ //for QUADQUADs and RQs.
 //                edgeInfo[2] = 0;    // means right side of edge is out & normal is OUTWARDS pointing.
 //            }
 //            edgeInfo[3] = 0 + j*(nX+1)+ nX*(nY+1); // general map
 //        edgeInfo[3] = 0; // local map
             edgeInfo[3] = -1; //means for periodic bc not working until set properly...
             if (edgeInfo[1]!=Q2I && edgeInfo[1]!= T2I && edgeInfo[1]!= Q2Q){
                 edgesAtBC.push_back(edgeInfo);
                 for (int icell = 0; icell < ngc; icell++){
                     std::vector<int> bm(3);
                     bcMap.push_back(nX[iii] + ngc + icell + (j + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     if (edgeInfo[1] == ABSORBINGBC){
                         bcMap.push_back(nX[iii] + ngc - 1 + (j + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     }else{
                         bcMap.push_back(nX[iii] + ngc - icell - 1 + (j + ngc)*(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
                     }
 
                     bcMap.push_back(edgesAtBC.size()-1);
                 }
             }else if (edgeInfo[1]==Q2Q){
                 for (int iii2 = iii+1; iii2 < nDomainsRQ + nDomainsQUADQUAD + nDomainsQUAD; iii2++){
                     if (iii2!=iii){
                         if (std::fabs(yMin[iii] - yMin[iii2])<eps0&&std::fabs(yMax[iii] - yMax[iii2])<eps0&&std::fabs(xMin[iii2] + 2*ngc*dX[iii] - xMax[iii])<eps0){
                             edgeInfo[3] = 0 + (j)*(nX[iii2] + 1) + nX[iii2]*(nY[iii2]+1) + submeshEdg0[iii2 + 1];
                             edgesAtBC.push_back(edgeInfo);
                             edgeInfo[3] = edgeInfo[0];
                             edgeInfo[0] = 0 + (j)*(nX[iii2] + 1) + nX[iii2]*(nY[iii2]+1) + submeshEdg0[iii2 + 1];
 //                            edgeInfo[2] = (edgeInfo[2] + 1)%2;
                             edgesAtBC.push_back(edgeInfo);
                             break;
                         }
                     }
                 }
             }else{
                 int nEdg_temp = edg.size();
                 for (int iEdg = 0; iEdg < nEdg_temp; iEdg++){
                     if(std::fabs((i + ngc)*dX[iii] + xMin[iii] - (xCoord[edg[iEdg][0]][0]+xCoord[edg[iEdg][1]][0])/2)<eps0 &&
                        std::fabs((j + 0.5 + ngc)*dY[iii] + yMin[iii] - (xCoord[edg[iEdg][0]][1]+xCoord[edg[iEdg][1]][1])/2)<eps0){
                         edgeInfo[3] = iEdg;
                         edgesAtBC.push_back(edgeInfo);
                         edgeInfo[3] = edgeInfo[0];
                         edgeInfo[0] = iEdg;
                         edgeInfo[2] = (edgeInfo[2] + 1)%2;
                         edgesAtBC.push_back(edgeInfo);
                         break;
                     }
                 }
             }
         }
     }
 
 //    for (int iiiii = 0 ; iiiii < edgesAtBC.size(); iiiii++){
 //
 //        if (edgesAtBC[iiiii][0] == 2007){
 //            std::cout << iiiii;
 //        }
 //    }
 
 
     // Create "right" cells of each edge and create edgeAtBCMap
     gmsh::logger::write("Create right cells of each edge and create edgeAtBCMap.");
     // Note: for each edge "left" cell is always in domain.
     nFaces = e2c.size();
     edgeInfo[3] = -1;
     int count = 0;
     for (int i = 0; i < nFaces; i++){
         if (e2c[i][1] == -1){
             edgeInfo[0] = i;
 //            if (edgeInfo[0] == 2007){
 //                std::cout << i;
 //            }
             int j = -1;
             count = 0;
             while (j < nBC-1&& count<2){
                 j++;
                 int n = listPhysicalNodeTagsEdges[j].size();
 
                 int ii = -1;
                 while (ii < n-1&& count < 2){
                     ii++;
                     if ((edg[i][0] == listPhysicalNodeTagsEdges[j][ii])|| (edg[i][1] == listPhysicalNodeTagsEdges[j][ii])){
                         count++;
                     }
                 }
                 if (count == 2){edgeInfo[1] = int(listPhysicalNameEdges[j]);}
             }
             edgeInfo[2] = 0;
             if (count == 2 && edgeInfo[1]!=Q2I && edgeInfo[1]!= T2I && edgeInfo[1]!= Q2Q){
                 edgesAtBC.push_back(edgeInfo);
             }
             count = 0;
         }
     }
 
 
     nCellsTot = submeshCell0.back();
 
     nedges = edg.size();
 
 
 
     switch (config.getModel()) {
         case EULER:
             nelem = 4;
             break;
         case BURGERS:
             nelem = 1;
             break;
             break;
         case SW:
             nelem = 3;
             break;
         case LINEARADVECTION:
             nelem = 1;
             break;
         case KPP:
             nelem = 1;
             break;
         case KPPY:
             nelem = 1;
             break;
         case KPPX:
             nelem = 1;
             break;
         default:
             assert(false && "FastHybridMesh: model not defined");
             break;
     };
 
     gmsh::logger::write("Caclulate normals of each edge.");
     normals.resize(submeshEdg0.back(),std::vector<double>(ndims));
     std::vector<double> tmpD(ndims);
     double m;
     for (int i = 0; i < nedges; i++){
 //        if (i == 14){
 //            std::cout << i << "\n";
 //        }
         tmpD[0] = xCoord[edg[i][1]][0]-xCoord[edg[i][0]][0];
         tmpD[1] = xCoord[edg[i][1]][1]-xCoord[edg[i][0]][1];
         m = std::sqrt(tmpD[0]*tmpD[0]+tmpD[1]*tmpD[1]);
         normals[i][0] = tmpD[1]/m;
         normals[i][1] = -tmpD[0]/m;
     }
     // calculate normals of edges for QUADQUADs, RQs and QUADs
     for (int ii = 0; ii < nDomainsQUADQUAD + nDomainsRQ + nDomainsQUAD; ii++){
         for (int i = 0; i < nX[ii]; i++){
             for (int j = 0; j < nY[ii]; j++){
                 normals[i +  j*(nX[ii]) + submeshEdg0[ii +1]][0] = 0;
                 normals[i +  j*(nX[ii]) + submeshEdg0[ii +1]][1] = 1;
             }
             int j = nY[ii];
             normals[i +  j*(nX[ii]) + submeshEdg0[ii +1]][0] = 0;
             normals[i +  j*(nX[ii]) + submeshEdg0[ii +1]][1] = -1;
         }
 
         for (int i = 0; i < nX[ii]; i++){
             for (int j = 0; j < nY[ii]; j++){
                 normals[i +  j*(nX[ii]+1) + nX[ii]*(nY[ii]+1) + submeshEdg0[ii +1]][0] = 1;
                 normals[i +  j*(nX[ii]+1) + nX[ii]*(nY[ii]+1) + submeshEdg0[ii +1]][1] = 0;
             }
         }
         int i = nX[ii];
         for (int j = 0; j < nY[ii]; j++){
             normals[i +  j*(nX[ii]+1) + nX[ii]*(nY[ii]+1) + submeshEdg0[ii +1]][0] = -1;
             normals[i +  j*(nX[ii]+1) + nX[ii]*(nY[ii]+1) + submeshEdg0[ii +1]][1] = 0;
         }
     }
 
 
     size.resize(nCellsTot);
     dS.resize(nCellsTot);
 
     double s;
     int nCellsTotTRIQUAD = submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT];
     std::vector<double> inRadius(nCellsTotTRIQUAD);
 //    std::cout << "test" << nCellsTot << " \n";
     double alpha = 0;
 //    int nNodes;
 
     for (int i = 0; i < nCellsTotTRIQUAD; i++){
         nNodes = cells[i].size();
         dS[i].resize(nNodes);
         if (nNodes == 3) {
             size[i] = 0.5 * std::fabs((xCoord[cells[i][1]][0] - xCoord[cells[i][0]][0]) *
                                       (xCoord[cells[i][2]][1] - xCoord[cells[i][0]][1])
                                       - (xCoord[cells[i][2]][0] - xCoord[cells[i][0]][0]) *
                                         (xCoord[cells[i][1]][1] - xCoord[cells[i][0]][1]));
 
             s = 0.5 * (sqrt(pow(xCoord[cells[i][1]][0] - xCoord[cells[i][0]][0], 2) +
                             pow(xCoord[cells[i][1]][1] - xCoord[cells[i][0]][1], 2))
                        + sqrt(pow(xCoord[cells[i][2]][0] - xCoord[cells[i][1]][0], 2) +
                               pow(xCoord[cells[i][2]][1] - xCoord[cells[i][1]][1], 2))
                        + sqrt(pow(xCoord[cells[i][0]][0] - xCoord[cells[i][2]][0], 2) +
                               pow(xCoord[cells[i][0]][1] - xCoord[cells[i][2]][1], 2)));
             inRadius[i] = size[i] / s;
             alpha = max(alpha, sqrt(pow(xCoord[cells[i][1]][0] - xCoord[cells[i][0]][0], 2) +
                                     pow(xCoord[cells[i][1]][1] - xCoord[cells[i][0]][1], 2)) / s);
             alpha = max(alpha, sqrt(pow(xCoord[cells[i][2]][0] - xCoord[cells[i][1]][0], 2) +
                                     pow(xCoord[cells[i][2]][1] - xCoord[cells[i][1]][1], 2)) / s);
             alpha = max(alpha, sqrt(pow(xCoord[cells[i][0]][0] - xCoord[cells[i][2]][0], 2) +
                                     pow(xCoord[cells[i][0]][1] - xCoord[cells[i][2]][1], 2)) / s);
         }else{
             size[i] = 0.5 * std::fabs((xCoord[cells[i][0]][0] - xCoord[cells[i][2]][0]) *
                                       (xCoord[cells[i][1]][1] - xCoord[cells[i][3]][1])
                                       - (xCoord[cells[i][0]][0] - xCoord[cells[i][2]][0]) *
                                         (xCoord[cells[i][3]][1] - xCoord[cells[i][1]][1]));
 
             s = 0.5 * (sqrt(pow(xCoord[cells[i][1]][0] - xCoord[cells[i][0]][0], 2) +
                             pow(xCoord[cells[i][1]][1] - xCoord[cells[i][0]][1], 2))
                        + sqrt(pow(xCoord[cells[i][2]][0] - xCoord[cells[i][1]][0], 2) +
                               pow(xCoord[cells[i][2]][1] - xCoord[cells[i][1]][1], 2))
                        + sqrt(pow(xCoord[cells[i][3]][0] - xCoord[cells[i][2]][0], 2) +
                               pow(xCoord[cells[i][3]][1] - xCoord[cells[i][2]][1], 2))
                        + sqrt(pow(xCoord[cells[i][0]][0] - xCoord[cells[i][3]][0], 2) +
                               pow(xCoord[cells[i][0]][1] - xCoord[cells[i][3]][1], 2)));
             inRadius[i] = size[i] / s;
             alpha = max(alpha, sqrt(pow(xCoord[cells[i][1]][0] - xCoord[cells[i][0]][0], 2) +
                                     pow(xCoord[cells[i][1]][1] - xCoord[cells[i][0]][1], 2)) / s);
             alpha = max(alpha, sqrt(pow(xCoord[cells[i][2]][0] - xCoord[cells[i][1]][0], 2) +
                                     pow(xCoord[cells[i][2]][1] - xCoord[cells[i][1]][1], 2)) / s);
             alpha = max(alpha, sqrt(pow(xCoord[cells[i][3]][0] - xCoord[cells[i][2]][0], 2) +
                                     pow(xCoord[cells[i][3]][1] - xCoord[cells[i][2]][1], 2)) / s);
             alpha = max(alpha, sqrt(pow(xCoord[cells[i][0]][0] - xCoord[cells[i][3]][0], 2) +
                                     pow(xCoord[cells[i][0]][1] - xCoord[cells[i][3]][1], 2)) / s);
         }
         for (int j = 0; j < nNodes; j++) {
             tmpD[0] = xCoord[cells[i][(j+1)%nNodes]][0] - xCoord[cells[i][j]][0];
             tmpD[1] = xCoord[cells[i][(j+1)%nNodes]][1] - xCoord[cells[i][j]][1];
             dS[i][j] = (tmpD[1] * normals[c2e[i][j]][0] - tmpD[0] * normals[c2e[i][j]][1]) / size[i];
         }
     }
     dxMin.resize(2);
     if (nCellsTotTRIQUAD>0) {
         double minInradius = *std::min_element(std::begin(inRadius), std::end(inRadius));
         dxMin[0] = minInradius / alpha;
         dxMin[1] = minInradius / alpha;
     }else{
         dxMin[0] = dX[0] / sqrt(2);
         dxMin[1] = dY[0] / sqrt(2);
     }
 
     for (int ii = 0; ii < nDomainsQUAD + nDomainsRQ + nDomainsQUADQUAD; ii++) {
         dxMin[0] = std::min(dxMin[0], dX[ii] / sqrt(2));
         dxMin[1] = std::min(dxMin[1], dY[ii] / sqrt(2));
     }
 
     std::vector<int> edg0(edg.size(), -1);
     int nedg;
 
 
 
     gmsh::logger::write("Set boundary conditions.");
     bc->initialize(this, config);
     gmsh::logger::write("Hybrid Mesh initialized.");
 }
 
 
 std::vector<double> FastHybridMesh::MidPoint(int i){
     std::vector<double> midPts(ndims, 0);
     if (i < cells.size()) {
         int nNodes = cells[i].size();
         for (int j = 0; j < nNodes; j++) {
             midPts[0] += xCoord[cells[i][j]][0];
             midPts[1] += xCoord[cells[i][j]][1];
         }
         midPts[0] /= nNodes;
         midPts[1] /= nNodes;
     }else{
         int i_test = 1;
         while(submeshCell0[i_test+1] > i){
             i_test++;
         }
         i -= submeshCell0[i_test];
         int i0 = i % (nX[i_test]+2*ngc);
         int j0 = (int) (i / (nX[i_test]+2*ngc));
         midPts[0] = ((double)i0 +0.5)*dX[i_test] + xMin[i_test];
         midPts[1] = ((double)j0 +0.5)*dY[i_test] + yMin[i_test];
     }
     return midPts;
 };
 
 void FastHybridMesh::netFlux(std::vector<std::vector<double>> &numFlux, NetFlux* numFluxFct, std::vector<std::vector<double>> &U,
                          std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, double t, double dt){
 
     // TRI, TRIQUAD, RQ part
     omp_set_num_threads(nThreads);
 #pragma omp parallel for
 
     for (int i = 0; i < nCells; i++){
         int ind_internal = internal[i];
 //        if (ind_internal == 1428){
 //            std::cout << i << "\n";
 //        }
 
 
         for (int k = 0; k < nelem; k++){
             numFlux[ind_internal][k] = 0;
         }
 
         int nNodes = c2e[ind_internal].size();
 
         for (int l = 0; l < nNodes; l++){
             int edge = c2e[ind_internal][l];
             std::vector<double> flux = numFluxFct->Feval(normals[edge], UReconstr[1][edge],UReconstr[0][edge],t,dt, dxMin);
             for (int k = 0; k < nelem; k++){
                 numFlux[ind_internal][k] -= flux[k]*dS[ind_internal][l];
             }
         }
     }
 
     // structured QUAD part
 //    std::vector<double> n_x{1,0};
 //    std::vector<double> n_y{0,1};
     for (int ii = 0; ii < nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ; ii++) {
 #pragma omp parallel for
         for (int i = 0; i < nX[ii]; i++) {
             for (int j = 0; j < nY[ii]; j++) {
                 for (int k = 0; k < nelem; k++) {
                     numFlux[i + ngc + (j + ngc) * (nX[ii] + 2 * ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k] = 0;
                 }
                 // calculate flux in y-direction
                 std::vector<double> flux1 = numFluxFct->Fevaly(normals[i + (j + 1) * nX[ii] + submeshEdg0[ii + 1]], UReconstr[1][i + (j + 1) * nX[ii] + submeshEdg0[ii + 1]],
                         UReconstr[0][i + (j + 1) * nX[ii] + submeshEdg0[ii + 1]], t, dt, dxMin);
                 std::vector<double> flux2 = numFluxFct->Fevaly(normals[i + (j) * nX[ii] + submeshEdg0[ii + 1]], UReconstr[1][i + (j) * nX[ii]+ submeshEdg0[ii + 1]],
                         UReconstr[0][i + (j) * nX[ii]+ submeshEdg0[ii + 1]],
                                                                t, dt, dxMin);
                 for (int k = 0; k < nelem; k++) {
                     numFlux[i + ngc + (j + ngc) * (nX[ii] + 2 * ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k] -= (flux1[k] - flux2[k]) / dY[ii];
                 }
                 //calculate flux in x-direction
                 flux1 = numFluxFct->Fevalx(normals[i + 1 + (j) * (nX[ii] + 1) + nX[ii] * (nY[ii] + 1)+ submeshEdg0[ii + 1]], UReconstr[1][i + 1 + (j) * (nX[ii] + 1) + nX[ii] * (nY[ii] + 1)+ submeshEdg0[ii + 1]],
                                            UReconstr[0][i + 1 + (j) * (nX[ii] + 1) + nX[ii] * (nY[ii] + 1)+ submeshEdg0[ii + 1]], t, dt, dxMin);
                 flux2 = numFluxFct->Fevalx(normals[i + (j) * (nX[ii] + 1) + nX[ii] * (nY[ii] + 1)+ submeshEdg0[ii + 1]], UReconstr[1][i + (j) * (nX[ii] + 1) + nX[ii] * (nY[ii] + 1)+ submeshEdg0[ii + 1]],
                                            UReconstr[0][i + (j) * (nX[ii] + 1) + nX[ii] * (nY[ii] + 1)+ submeshEdg0[ii + 1]], t, dt, dxMin);
                 for (int k = 0; k < nelem; k++) {
                     numFlux[i + ngc + (j + ngc) * (nX[ii] + 2 * ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k] -= (flux1[k] - flux2[k]) / dX[ii];
                 }
             }
         }
     }
 };
 
 std::vector<std::vector<double> > FastHybridMesh::evalFunc(AbstractFunct *Funct, evaluationMethod evalMethod){
 
     std::vector<std::vector<double> > ret(nCellsTot, std::vector<double>(nelem,0));
     switch (evalMethod){
         case MEANVALUE:{
 
             std::vector<double> node_xy2;
             int i;
             int order_num;
             int rule;
             double *w;
             double *xy_dunavant;
 //            double *xy2;
 
             std::vector<double> weights_legendre;
             std::vector<double> x_legendre;
             std::vector<double> jacob;
 
             rule = 2;
 
             int n = rule;
             x_legendre.resize(n);
             weights_legendre.resize(n);
             double *w0 = new double[n];
             double *x = new double[n];
 
             legendre_set(n,x,w0);
             if (n < 2){
                 x_legendre[0] = 0;
                 weights_legendre[0] = 2;
             }else {
                 for (int i = 0; i < n; i++) {
                     x_legendre[i] = x[i];
                     weights_legendre[i] = w0[i];
                 }
             }
 
 
             order_num = dunavant_order_num ( rule );
             xy_dunavant = new double[2*order_num];
             std::vector<double> xy2;
             w = new double[order_num];
 
             dunavant_rule ( rule, order_num, xy_dunavant, w );
             std::vector<double> x_ev(2);
             std::vector<double> F_ev;
             std::vector<int> cell_i;
             int nNodes;
             int nCellsTotTRIQUAD = submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT];
 
             for (i = 0; i < nCellsTotTRIQUAD; i++){
                 cell_i = cells[i];
                 nNodes = cell_i.size();
                 node_xy2.resize(2*nNodes);
 
                 for (int j = 0; j < nNodes; j++) {
                     node_xy2[2*j] = xCoord[cells[i][j]][0];
                     node_xy2[2*j + 1] = xCoord[cells[i][j]][1];
                 }
 
                 if (nNodes == 3) {
                     xy2.resize(2*order_num);
                     jacob.resize(order_num, 1.0);
                     std::fill(jacob.begin(), jacob.end(), 1.0);
 
                     reference_to_physical_t3 ( node_xy2, order_num, xy_dunavant, xy2 );
                     for (int j = 0; j < order_num; j++){
                         jacob[j] *= w[j];
                     }
 //                    std::fill(ret[i].begin(), ret[i].end(), 0);
 //                    for (int k = 0; k<order_num; k++){
 //                        x_ev[0] = xy2[0+k*2];
 //                        x_ev[1] = xy2[1+k*2];
 //
 //                        F_ev = Funct->F(x_ev);
 //                        for (int l = 0; l<nelem; l++){
 //                            ret[i][l] += w[k]*F_ev[l];
 //                        }
 //                    }
                 } else{
 //            assert(false && "nNodes == 4, not implemented, yet!");
                     xy2.resize(2*n*n);
                     jacob.resize(n*n);
                     reference_to_physical_t4(node_xy2, n, x_legendre, weights_legendre, xy2, jacob);
                 }
 
 
                 int nI = jacob.size();
                 std::fill(ret[i].begin(), ret[i].end(), 0);
                 for (int k = 0; k<nI; k++){
                     x_ev[0] = xy2[0+k*2];
                     x_ev[1] = xy2[1+k*2];
                     F_ev = Funct->F(x_ev);
                     for (int l = 0; l<nelem; l++){
                         ret[i][l] += jacob[k]*F_ev[l];
                     }
                 }
             }
 
             // evaluate QUAD part
 //            rule = 4;
 
             n = rule;
             x_legendre.resize(n);
             weights_legendre.resize(n);
 //            double *w0 = new double[n];
 //            double *x = new double[n];
 
             legendre_set(n,x,w0);
             if (n < 2){
                 x_legendre[0] = 0;
                 weights_legendre[0] = 2;
             }else {
                 for (int i = 0; i < n; i++) {
                     x_legendre[i] = x[i];
                     weights_legendre[i] = w0[i];
                 }
             }
 
 
             order_num = dunavant_order_num ( rule );
             xy_dunavant = new double[2*order_num];
 //            std::vector<double> xy2;
             w = new double[order_num];
 
             dunavant_rule ( rule, order_num, xy_dunavant, w );
 //            int nX = xCoord[0].size()-1;
 //            int nY = xCoord[1].size()-1;
 
             for (int ii = 0; ii < nDomainsQUAD + nDomainsRQ + nDomainsQUADQUAD; ii++) {
                 for (int i = 0; i < nX[ii] + 2 * ngc; i++) {
                     for (int j = 0; j < nY[ii] + 2 * ngc; j++) {
                         nNodes = 4;
                         node_xy2.resize(2 * nNodes);
 
                         node_xy2[2 * 0] = i*dX[ii] + xMin[ii];
                         node_xy2[2 * 0 + 1] = j*dY[ii] + yMin[ii];
                         node_xy2[2 * 1] = (i+1)*dX[ii] + xMin[ii];
                         node_xy2[2 * 1 + 1] = j*dY[ii] + yMin[ii];
                         node_xy2[2 * 2] = (i+1)*dX[ii] + xMin[ii];
                         node_xy2[2 * 2 + 1] = (j+1)*dY[ii] + yMin[ii];
                         node_xy2[2 * 3] = i*dX[ii] + xMin[ii];
                         node_xy2[2 * 3 + 1] = (j+1)*dY[ii] + yMin[ii];
 
                         xy2.resize(2 * n * n);
                         jacob.resize(n * n);
                         reference_to_physical_t4(node_xy2, n, x_legendre, weights_legendre, xy2, jacob);
 
                         int nI = jacob.size();
                         std::fill(ret[i + j * (nX[ii] + 2 * ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]].begin(),
                                 ret[i + j * (nX[ii] + 2 * ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]].end(), 0);
                         for (int k = 0; k < nI; k++) {
                             x_ev[0] = xy2[0 + k * 2];
                             x_ev[1] = xy2[1 + k * 2];
                             F_ev = Funct->F(x_ev);
                             for (int l = 0; l < nelem; l++) {
                                 ret[i + j * (nX[ii] + 2 * ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][l] += jacob[k] * F_ev[l];
                             }
                         }
                     }
                 }
             }
             delete [] w;
             delete [] xy_dunavant;
             delete [] x;
             delete [] w0;
             break;
         }
         case PTWISE:{
             for (int i = 0; i < nCellsTot; i++){
                 std::vector<double> m = MidPoint(i);
                 ret[i] = Funct->F(MidPoint(i));
             }
             break;
         }
 
     }
     return ret;
 };
 
 std::vector<std::vector<double>> FastHybridMesh::edge2Vertex(int i){
     std::vector<std::vector<double>> ret(2,std::vector<double>(2,0));
     std::vector<int> edg_i = edg[i];
     ret[0][0] = xCoord[edg_i[0]][0];
     ret[0][1] = xCoord[edg_i[1]][0];
     ret[1][0] = xCoord[edg_i[0]][1];
     ret[1][1] = xCoord[edg_i[1]][1];
     return ret;
 }
 
 std::vector<int> FastHybridMesh::getSubmeshCell0(){return submeshCell0;};
 
 std::vector<int> FastHybridMesh::getSubmeshEdg0(){return submeshEdg0;};
 
 std::vector<int> FastHybridMesh::getNx(){return nX;};
 
 std::vector<int> FastHybridMesh::getNy(){return nY;};
 
 std::vector<double> FastHybridMesh::getDx(){return dX;};
 
 std::vector<double> FastHybridMesh::getDy(){return dY;};
 
 std::vector<double> FastHybridMesh::getXmin(){return xMin;};
 
 std::vector<double> FastHybridMesh::getYmin(){return yMin;};
 
 void FastHybridMesh::convertQUAD2STRUCTURED(std::vector<int> &TRIQUAD2QUAD, double &eps0,
         std::vector<int> &coordmap, std::vector<int> &listCoordMap){
     xMin.resize(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ);
     xMax.resize(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ);
     yMin.resize(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ);
     yMax.resize(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ);
     dX.resize(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ);
     dY.resize(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ);
     nX.resize(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ);
     nY.resize(nDomainsQUAD + nDomainsQUADQUAD + nDomainsRQ);
 
     int iElement;
     for ( int i = 0; i < nDomainsQUAD; i++){
         // find xMin, xMax, yMin, yMax
         xMin[i] = xCoord[cells[submeshCell0[nDomainsTRI + i]][0]][0];
         xMax[i] = xCoord[cells[submeshCell0[nDomainsTRI + i]][0]][0];
         yMin[i] = xCoord[cells[submeshCell0[nDomainsTRI + i]][0]][1];
         yMax[i] = xCoord[cells[submeshCell0[nDomainsTRI + i]][0]][1];
 
         for (int ii = submeshCell0[nDomainsTRI + i]; ii < submeshCell0[nDomainsTRI + i + 1]; ii++){
 
             for (int j = 0; j < 4; j++){
                 if (xMin[i] >  xCoord[cells[ii][j]][0]){
                     xMin[i] = xCoord[cells[ii][j]][0];
                 }
                 if (xMax[i] <  xCoord[cells[ii][j]][0]){
                     xMax[i] = xCoord[cells[ii][j]][0];
                 }
                 if (yMin[i] >  xCoord[cells[ii][j]][1]){
                     yMin[i] = xCoord[cells[ii][j]][1];
                 }
                 if (yMax[i] <  xCoord[cells[ii][j]][1]){
                     yMax[i] = xCoord[cells[ii][j]][1];
                 }
             }
         }
 
         iElement = submeshCell0[nDomainsTRI + i];
         dX[i] =  std:: max(std::fabs(xCoord[cells[iElement][0]][0] - xCoord[cells[iElement][1]][0])
                 , std::fabs(xCoord[cells[iElement][1]][0] - xCoord[cells[iElement][2]][0]));
         dY[i] =  std:: max(std::fabs(xCoord[cells[iElement][0]][1] - xCoord[cells[iElement][1]][1])
                 , std::fabs(xCoord[cells[iElement][1]][1] - xCoord[cells[iElement][2]][1]));
         nX[i] = round((xMax[i]- xMin[i])/dX[i]);
         nY[i] = round((yMax[i]-yMin[i])/dY[i]);
 
         dX[i] = (xMax[i] - xMin[i])/(nX[i]);
         dY[i] = (yMax[i] - yMin[i])/(nY[i]);
         xMax[i] += ngc*dX[i];
         xMin[i] -= ngc*dX[i];
         yMax[i] += ngc*dY[i];
         yMin[i] -= ngc*dY[i];
 
 //        eps0 = std::min(std::min(dX[i]/10, dY[i]/10), eps0);
 //        coordMap[i].resize((nX[i] + 1)*(nY[i] + 1));
 //        listCoordMap[i + 1] = (listCoordMap.back() + (nX[i] + 1)*(nY[i] + 1));
 //        int nCoord = xCoord.size();
 //        for (int ii = 0; ii <= nX[i]; ii++){
 //            for (int jj = 0; jj <= nY[i]; jj++){
 //                for (int k = 0; k < nCoord; k++){
 //                    if ((std::abs(xCoord[k][0]- (ii*dX[i] + xMin[i]+ngc*dX[i]))< eps0) &&
 //                        (std::abs(xCoord[k][1]- (jj*dY[i] + yMin[i]+ngc*dY[i]))< eps0)){
 //                        coordMap[i][ii + (jj )*(nX[i] + 1 )] = k ;
 //                        coordmap[listCoordMap[i] + ii + (jj )*(nX[i] + 1 )] = k;
 //                    }
 //                }
 //            }
 //        }
     }
 
     int i0;
     for ( int i = 0; i < nDomainsQUADQUAD; i++){
         // find xMin, xMax, yMin, yMax
         i0 = nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT;
         xMin[i + nDomainsQUAD] = xCoord[cells[submeshCell0[i0 + i]][0]][0];
         xMax[i + nDomainsQUAD] = xCoord[cells[submeshCell0[i0 + i]][0]][0];
         yMin[i + nDomainsQUAD] = xCoord[cells[submeshCell0[i0 + i]][0]][1];
         yMax[i + nDomainsQUAD] = xCoord[cells[submeshCell0[i0 + i]][0]][1];
 
         for (int ii = submeshCell0[i0 + i]; ii < submeshCell0[i0 + i + 1]; ii++){
 
             for (int j = 0; j < 4; j++){
                 if (xMin[i + nDomainsQUAD] >  xCoord[cells[ii][j]][0]){
                     xMin[i + nDomainsQUAD] = xCoord[cells[ii][j]][0];
                 }
                 if (xMax[i + nDomainsQUAD] <  xCoord[cells[ii][j]][0]){
                     xMax[i + nDomainsQUAD] = xCoord[cells[ii][j]][0];
                 }
                 if (yMin[i + nDomainsQUAD] >  xCoord[cells[ii][j]][1]){
                     yMin[i + nDomainsQUAD] = xCoord[cells[ii][j]][1];
                 }
                 if (yMax[i + nDomainsQUAD] <  xCoord[cells[ii][j]][1]){
                     yMax[i + nDomainsQUAD] = xCoord[cells[ii][j]][1];
                 }
             }
         }
         iElement = submeshCell0[i0 + i];
         dX[i + nDomainsQUAD] =  std:: max(std::fabs(xCoord[cells[iElement][0]][0] - xCoord[cells[iElement][1]][0])
                 , std::fabs(xCoord[cells[iElement][1]][0] - xCoord[cells[iElement][2]][0]));
         dY[i + nDomainsQUAD] =  std:: max(std::fabs(xCoord[cells[iElement][0]][1] - xCoord[cells[iElement][1]][1])
                 , std::fabs(xCoord[cells[iElement][1]][1] - xCoord[cells[iElement][2]][1]));
         nX[i + nDomainsQUAD] = round((xMax[i + nDomainsQUAD]- xMin[i + nDomainsQUAD])/dX[i + nDomainsQUAD]);
         nY[i + nDomainsQUAD] = round((yMax[i + nDomainsQUAD]-yMin[i + nDomainsQUAD])/dY[i + nDomainsQUAD]);
 
         dX[i + nDomainsQUAD] = (xMax[i + nDomainsQUAD] - xMin[i + nDomainsQUAD])/(nX[i + nDomainsQUAD]);
         dY[i + nDomainsQUAD] = (yMax[i + nDomainsQUAD] - yMin[i + nDomainsQUAD])/(nY[i + nDomainsQUAD]);
         xMax[i + nDomainsQUAD] += ngc*dX[i + nDomainsQUAD];
         xMin[i + nDomainsQUAD] -= ngc*dX[i + nDomainsQUAD];
         yMax[i + nDomainsQUAD] += ngc*dY[i + nDomainsQUAD];
         yMin[i + nDomainsQUAD] -= ngc*dY[i + nDomainsQUAD];
 
 //        eps0 = std::min(std::min(dX[i + nDomainsQUAD]/10, dY[i + nDomainsQUAD]/10), eps0);
 //        coordMap[i + nDomainsQUAD].resize((nX[i + nDomainsQUAD] + 1)*(nY[i + nDomainsQUAD] + 1));
 //        listCoordMap[i + nDomainsQUAD + 1] = (nX[i + nDomainsQUAD] + 1)*(nY[i + nDomainsQUAD] + 1);
 //        int nCoord = xCoord.size();
 //        for (int ii = 0; ii <= nX[i + nDomainsQUAD]; ii++){
 //            for (int jj = 0; jj <= nY[i + nDomainsQUAD]; jj++){
 //                for (int k = 0; k < nCoord; k++){
 //                    if ((std::abs(xCoord[k][0]- (ii*dX[i + nDomainsQUAD] + xMin[i + nDomainsQUAD]+ngc*dX[i + nDomainsQUAD]))< eps0) &&
 //                        (std::abs(xCoord[k][1]- (jj*dY[i + nDomainsQUAD] + yMin[i + nDomainsQUAD]+ngc*dY[i + nDomainsQUAD]))< eps0)){
 //                        coordMap[i + nDomainsQUAD][ii + (jj )*(nX[i + nDomainsQUAD] + 1 )] = k ;
 //                        coordmap[listCoordMap[i + nDomainsQUAD] + ii + (jj )*(nX[i + nDomainsQUAD] + 1 )] = k;
 //                    }
 //                }
 //            }
 //        }
     }
 
     for ( int i = 0; i < nDomainsRQ; i++){
         // find xMin, xMax, yMin, yMax
         i0 = nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT + nDomainsQUADQUAD;
         xMin[i + nDomainsQUAD + nDomainsQUADQUAD] = xCoord[cells[submeshCell0[i0 + i]][0]][0];
         xMax[i + nDomainsQUAD + nDomainsQUADQUAD] = xCoord[cells[submeshCell0[i0 + i]][0]][0];
         yMin[i + nDomainsQUAD + nDomainsQUADQUAD] = xCoord[cells[submeshCell0[i0 + i]][0]][1];
         yMax[i + nDomainsQUAD + nDomainsQUADQUAD] = xCoord[cells[submeshCell0[i0 + i]][0]][1];
 
         for (int ii = submeshCell0[i0 + i]; ii < submeshCell0[i0 + i + 1]; ii++){
 
             for (int j = 0; j < 4; j++){
                 if (xMin[i + nDomainsQUAD + nDomainsQUADQUAD] >  xCoord[cells[ii][j]][0]){
                     xMin[i + nDomainsQUAD + nDomainsQUADQUAD] = xCoord[cells[ii][j]][0];
                 }
                 if (xMax[i + nDomainsQUAD + nDomainsQUADQUAD] <  xCoord[cells[ii][j]][0]){
                     xMax[i + nDomainsQUAD + nDomainsQUADQUAD] = xCoord[cells[ii][j]][0];
                 }
                 if (yMin[i + nDomainsQUAD + nDomainsQUADQUAD] >  xCoord[cells[ii][j]][1]){
                     yMin[i + nDomainsQUAD + nDomainsQUADQUAD] = xCoord[cells[ii][j]][1];
                 }
                 if (yMax[i + nDomainsQUAD + nDomainsQUADQUAD] <  xCoord[cells[ii][j]][1]){
                     yMax[i + nDomainsQUAD + nDomainsQUADQUAD] = xCoord[cells[ii][j]][1];
                 }
             }
         }
         iElement = submeshCell0[i0 + i];
         dX[i + nDomainsQUAD + nDomainsQUADQUAD] =  std:: max(std::fabs(xCoord[cells[iElement][0]][0] - xCoord[cells[iElement][1]][0])
                 , std::fabs(xCoord[cells[iElement][1]][0] - xCoord[cells[iElement][2]][0]));
         dY[i + nDomainsQUAD + nDomainsQUADQUAD] =  std:: max(std::fabs(xCoord[cells[iElement][0]][1] - xCoord[cells[iElement][1]][1])
                 , std::fabs(xCoord[cells[iElement][1]][1] - xCoord[cells[iElement][2]][1]));
         nX[i + nDomainsQUAD + nDomainsQUADQUAD] = round((xMax[i + nDomainsQUAD + nDomainsQUADQUAD]- xMin[i + nDomainsQUAD + nDomainsQUADQUAD])/dX[i + nDomainsQUAD + nDomainsQUADQUAD]);
         nY[i + nDomainsQUAD + nDomainsQUADQUAD] = round((yMax[i + nDomainsQUAD + nDomainsQUADQUAD]-yMin[i + nDomainsQUAD + nDomainsQUADQUAD])/dY[i + nDomainsQUAD + nDomainsQUADQUAD]);
 
         dX[i + nDomainsQUAD + nDomainsQUADQUAD] = (xMax[i + nDomainsQUAD + nDomainsQUADQUAD] - xMin[i + nDomainsQUAD + nDomainsQUADQUAD])/(nX[i + nDomainsQUAD + nDomainsQUADQUAD]);
         dY[i + nDomainsQUAD + nDomainsQUADQUAD] = (yMax[i + nDomainsQUAD + nDomainsQUADQUAD] - yMin[i + nDomainsQUAD + nDomainsQUADQUAD])/(nY[i + nDomainsQUAD + nDomainsQUADQUAD]);
         xMax[i + nDomainsQUAD + nDomainsQUADQUAD] += ngc*dX[i + nDomainsQUAD + nDomainsQUADQUAD];
         xMin[i + nDomainsQUAD + nDomainsQUADQUAD] -= ngc*dX[i + nDomainsQUAD + nDomainsQUADQUAD];
         yMax[i + nDomainsQUAD + nDomainsQUADQUAD] += ngc*dY[i + nDomainsQUAD + nDomainsQUADQUAD];
         yMin[i + nDomainsQUAD + nDomainsQUADQUAD] -= ngc*dY[i + nDomainsQUAD + nDomainsQUADQUAD];
 
 //        eps0 = std::min(std::min(dX[i + nDomainsQUAD + nDomainsQUADQUAD]/10, dY[i + nDomainsQUAD + nDomainsQUADQUAD]/10), eps0);
 //        coordMap[i + nDomainsQUAD + nDomainsQUADQUAD].resize((nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 1)*(nY[i + nDomainsQUAD + nDomainsQUADQUAD] + 1));
 //        listCoordMap[i + nDomainsQUAD + nDomainsQUADQUAD + 1] = (nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 1)*(nY[i + nDomainsQUAD + nDomainsQUADQUAD] + 1);
 //        int nCoord = xCoord.size();
 //        for (int ii = 0; ii <= nX[i + nDomainsQUAD + nDomainsQUADQUAD]; ii++){
 //            for (int jj = 0; jj <= nY[i + nDomainsQUAD + nDomainsQUADQUAD]; jj++){
 //                for (int k = 0; k < nCoord; k++){
 //                    if ((std::abs(xCoord[k][0]- (ii*dX[i + nDomainsQUAD + nDomainsQUADQUAD] + xMin[i + nDomainsQUAD + nDomainsQUADQUAD]+ngc*dX[i + nDomainsQUAD + nDomainsQUADQUAD]))< eps0) &&
 //                        (std::abs(xCoord[k][1]- (jj*dY[i + nDomainsQUAD + nDomainsQUADQUAD] + yMin[i + nDomainsQUAD + nDomainsQUADQUAD]+ngc*dY[i + nDomainsQUAD + nDomainsQUADQUAD]))< eps0)){
 //                        coordMap[i + nDomainsQUAD + nDomainsQUADQUAD][ii + (jj )*(nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 1 )] = k ;
 //                        coordmap[listCoordMap[i + nDomainsQUAD + nDomainsQUADQUAD] + ii + (jj )*(nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 1 )] = k;
 //                    }
 //                }
 //            }
 //        }
     }
 
     int NN = nX.size();
     int nEdges = 0;
     for (int i = 0; i< NN; i++){
         nEdges += (nX[i]+1)*(nY[i]+1);
     }
     coordmap.resize(nEdges);
 
     for ( int i = 0; i < nDomainsQUADQUAD; i++) {
         eps0 = std::min(std::min(dX[i + nDomainsQUAD]/10, dY[i + nDomainsQUAD]/10), eps0);
 //        coordMap[i + nDomainsQUAD].resize((nX[i + nDomainsQUAD] + 1)*(nY[i + nDomainsQUAD] + 1));
         listCoordMap[i + 1] = listCoordMap[i] + (nX[i + nDomainsQUAD] + 1)*(nY[i + nDomainsQUAD] + 1);
         int nCoord = xCoord.size();
         for (int ii = 0; ii <= nX[i + nDomainsQUAD]; ii++){
             for (int jj = 0; jj <= nY[i + nDomainsQUAD]; jj++){
                 for (int k = 0; k < nCoord; k++){
                     if ((std::abs(xCoord[k][0]- (ii*dX[i + nDomainsQUAD] + xMin[i + nDomainsQUAD]+ngc*dX[i + nDomainsQUAD]))< eps0) &&
                         (std::abs(xCoord[k][1]- (jj*dY[i + nDomainsQUAD] + yMin[i + nDomainsQUAD]+ngc*dY[i + nDomainsQUAD]))< eps0)){
 //                        coordMap[i + nDomainsQUAD][ii + (jj )*(nX[i + nDomainsQUAD] + 1 )] = k ;
                         coordmap[listCoordMap[i] + ii + (jj )*(nX[i + nDomainsQUAD] + 1 )] = k;
                     }
                 }
             }
         }
     }
 
     for ( int i = 0; i < nDomainsRQ; i++){
         eps0 = std::min(std::min(dX[i + nDomainsQUAD + nDomainsQUADQUAD]/10, dY[i + nDomainsQUAD + nDomainsQUADQUAD]/10), eps0);
 //        coordMap[i + nDomainsQUAD + nDomainsQUADQUAD].resize((nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 1)*(nY[i + nDomainsQUAD + nDomainsQUADQUAD] + 1));
         listCoordMap[i + nDomainsQUADQUAD + 1] = listCoordMap[i + nDomainsQUADQUAD] + (nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 1)*(nY[i + nDomainsQUAD + nDomainsQUADQUAD] + 1);
         int nCoord = xCoord.size();
         for (int ii = 0; ii <= nX[i + nDomainsQUAD + nDomainsQUADQUAD]; ii++){
             for (int jj = 0; jj <= nY[i + nDomainsQUAD + nDomainsQUADQUAD]; jj++){
                 for (int k = 0; k < nCoord; k++){
                     if ((std::abs(xCoord[k][0]- (ii*dX[i + nDomainsQUAD + nDomainsQUADQUAD] + xMin[i + nDomainsQUAD + nDomainsQUADQUAD]+ngc*dX[i + nDomainsQUAD + nDomainsQUADQUAD]))< eps0) &&
                         (std::abs(xCoord[k][1]- (jj*dY[i + nDomainsQUAD + nDomainsQUADQUAD] + yMin[i + nDomainsQUAD + nDomainsQUADQUAD]+ngc*dY[i + nDomainsQUAD + nDomainsQUADQUAD]))< eps0)){
 //                        coordMap[i + nDomainsQUAD + nDomainsQUADQUAD][ii + (jj )*(nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 1 )] = k ;
                         coordmap[listCoordMap[i + nDomainsQUADQUAD] + ii + (jj )*(nX[i + nDomainsQUAD + nDomainsQUADQUAD] + 1 )] = k;
                     }
                 }
             }
         }
     }
 
     for ( int i = 0; i < nDomainsQUAD; i++){
         eps0 = std::min(std::min(dX[i]/10, dY[i]/10), eps0);
 //        coordMap[i].resize((nX[i] + 1)*(nY[i] + 1));
         listCoordMap[i + nDomainsRQ + nDomainsQUADQUAD + 1] = listCoordMap[i + nDomainsRQ + nDomainsQUADQUAD] + (listCoordMap.back() + (nX[i] + 1)*(nY[i] + 1));
         int nCoord = xCoord.size();
         for (int ii = 0; ii <= nX[i]; ii++){
             for (int jj = 0; jj <= nY[i]; jj++){
                 for (int k = 0; k < nCoord; k++){
                     if ((std::abs(xCoord[k][0]- (ii*dX[i] + xMin[i]+ngc*dX[i]))< eps0) &&
                         (std::abs(xCoord[k][1]- (jj*dY[i] + yMin[i]+ngc*dY[i]))< eps0)){
 //                        coordMap[i][ii + (jj )*(nX[i] + 1 )] = k ;
                         coordmap[listCoordMap[i + nDomainsRQ + nDomainsQUADQUAD] + ii + (jj )*(nX[i] + 1 )] = k;
                     }
                 }
             }
         }
     }
 }
 
 void FastHybridMesh::addGhostCells2QUAD(std::vector<MeshType> &listPhysicalNameElements, std::vector<std::vector<std::vector<int>>> &cells_add,
         double & eps0, std::vector<int> &TRIQUAD2QUAD, std::vector<std::vector<int>> &bcMap_loc) {
     // for each TRIQUAD
     for (int i = nDomainsQUAD + nDomainsTRI; i < nDomainsQUAD + nDomainsTRI + nDomainsTRIQUAD ;i++){
         double xMin_temp = xCoord[cells[submeshCell0[i]][0]][0];
         double xMax_temp = xCoord[cells[submeshCell0[i]][0]][0];
         double yMin_temp = xCoord[cells[submeshCell0[i]][0]][1];
         double yMax_temp = xCoord[cells[submeshCell0[i]][0]][1];
 
         // find xMin, xMax, yMin, yMax for TRIQUAD
         for (int ii = submeshCell0[i]; ii < submeshCell0[i + 1]; ii++){
 
 //            double eps0 = std::min(dX[i]/10, dY[i]/10);
 
             for (int j = 0; j < 4; j++){
                 if (xMin_temp >  xCoord[cells[ii][j]][0]){
                     xMin_temp = xCoord[cells[ii][j]][0];
                 }
                 if (xMax_temp <  xCoord[cells[ii][j]][0]){
                     xMax_temp = xCoord[cells[ii][j]][0];
                 }
                 if (yMin_temp >  xCoord[cells[ii][j]][1]){
                     yMin_temp = xCoord[cells[ii][j]][1];
                 }
                 if (yMax_temp <  xCoord[cells[ii][j]][1]){
                     yMax_temp = xCoord[cells[ii][j]][1];
                 }
             }
         }
         int nCellsTRIQUAD = submeshCell0[i+1]-submeshCell0[i];
         cells_add[i - (nDomainsQUAD + nDomainsTRI)].resize(nCellsTRIQUAD, std::vector<int>(4));
         std::vector<double> nodes0(8);
 
 
         double xM, yM, xM0, yM0;
         for (int ii = 0; ii < nDomainsQUAD; ii++) {
             if ((std::fabs(xMax_temp-xMax[ii])<eps0)&&(std::fabs(yMax_temp-yMax[ii]+ngc*dY[ii])<eps0)&&
             (std::fabs(yMin_temp-yMin[ii]-ngc*dY[ii])<eps0)){ //TRIQUAD domain i is right of QUAD domain ii
                 // enlarge QUAD domain ii and find Map from TRIQUAD to ghost cells of QUAD
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 TRIQUAD2QUAD[i - (nDomainsQUAD + nDomainsTRI)] = ii;
 
                 std::vector<int> bm(0);
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){
                     xM = 0;
                     yM = 0;
 //                    bm[1] = jj - submeshCell0[i];
                     for (int j = 0; j < 4; j++){
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int i0 = (int) std::round((xM - xMin[ii])/dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii])/dY[ii] - 0.5);
 
                     bm.push_back(i0 + j0*(2*ngc + nX[ii]));
                     bm.push_back(jj-submeshCell0[i]);
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
 
 
                 // enlarge TRIQUAD domain i
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){
 
                     eps0 = std::min(dX[ii]/10, dY[ii]/10);
 
                     xM = 0;
                     yM = 0;
                     for (int j = 0; j < 4; j++){
                         nodes0[2*j] = xCoord[cells[jj][j]][0]-(xMax_temp-xMin_temp);
                         nodes0[2*j + 1] = xCoord[cells[jj][j]][1];
                         xM += xCoord[cells[jj][j]][0]-(xMax_temp-xMin_temp);
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int cellFound = 0;
 
                     int l = submeshCell0[ii + nDomainsTRI];
                     bool foundCell = false;
                     while ((l < submeshCell0[ii + 1 + nDomainsTRI])&&((!foundCell)) ) {
                         int j = 0;
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++){
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
                         if ((std::abs(xM0 - xM) < eps0)&&(std::abs(yM0 - yM) < eps0)){
                             cells_add[i - (nDomainsQUAD + nDomainsTRI)][jj-submeshCell0[i]]= cells[l];
                             cellFound++;
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii])/dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii])/dY[ii] - 0.5);
                     // note: submeshCell0[ii + nDomainsTRI + nDomainsQUAD+1] - submeshCell0[ii + nDomainsTRI + nDomainsQUAD]
                     // is the number of cells in the current TRIQUAD ii.
                     bm.push_back(jj-submeshCell0[i] + submeshCell0[i + 1] - submeshCell0[i]);
                     bm.push_back(i0 + j0*(2*ngc + nX[ii]));
                     bm.push_back(-1); // note: -1 stands for interior update
 //                    assert(cellFound == 1 && "FastHybridMesh: Wrong number of cells found!");
                 }
                 bcMap_loc.push_back(bm);
                 break;
             }else if((std::fabs(xMin_temp - xMin[ii])<eps0)&&(std::fabs(yMax_temp-yMax[ii]+ngc*dY[ii])<eps0)&&
                      (std::fabs(yMin_temp-yMin[ii]-ngc*dY[ii])<eps0)){ //TRIQUAD domain i is left of QUAD domain ii
                 // enlarge QUAD domain ii
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 TRIQUAD2QUAD[i - (nDomainsQUAD + nDomainsTRI)] = ii;
 
                 std::vector<int> bm(0);
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){
                     xM = 0;
                     yM = 0;
 //                    bm[1] = jj - submeshCell0[i];
                     for (int j = 0; j < 4; j++){
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int i0 = (int) std::round((xM - xMin[ii])/dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii])/dY[ii] - 0.5);
 
                     bm.push_back(i0 + j0*(2*ngc + nX[ii]));
                     bm.push_back(jj-submeshCell0[i]);
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
 
                 // enlarge TRIQUAD domain i
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){
 
                     double eps0 = std::min(dX[ii]/10, dY[ii]/10);
                     xM = 0;
                     yM = 0;
                     for (int j = 0; j < 4; j++){
                         nodes0[2*j] = xCoord[cells[jj][j]][0]+(xMax_temp-xMin_temp);
                         nodes0[2*j + 1] = xCoord[cells[jj][j]][1];
                         xM += xCoord[cells[jj][j]][0]+(xMax_temp-xMin_temp);
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int cellFound = 0;
 
                     int l = submeshCell0[ii + nDomainsTRI];
                     bool foundCell = false;
                     while ((l < submeshCell0[ii + 1 + nDomainsTRI])&&((!foundCell)) ) {
                         int j = 0;
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++){
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
                         if ((std::abs(xM0 - xM) < eps0)&&(std::abs(yM0 - yM) < eps0)){
                             cells_add[i - (nDomainsQUAD + nDomainsTRI)][jj-submeshCell0[i]]= cells[l];
                             cellFound++;
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii])/dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii])/dY[ii] - 0.5);
 
                     bm.push_back(jj-submeshCell0[i] + submeshCell0[i + 1] - submeshCell0[i]);
                     bm.push_back(i0 + j0*(2*ngc + nX[ii]));
                     bm.push_back(-1); // note: -1 stands for interior update
 //                    assert(cellFound == 1 && "FastHybridMesh: Wrong number of cells found!");
                 }
                 bcMap_loc.push_back(bm);
                 break;
             }else if((std::fabs(yMax_temp - yMax[ii])<eps0)&&(std::fabs(xMax_temp-xMax[ii]+ngc*dX[ii])<eps0)&&
                      (std::fabs(xMin_temp-xMin[ii]-ngc*dX[ii])<eps0)){ //TRIQUAD domain i is on top of QUAD domain ii
                 // enlarge QUAD domain ii
                 assert((int) std::round((yMax_temp - yMin_temp)/dY[ii] )==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 TRIQUAD2QUAD[i - (nDomainsQUAD + nDomainsTRI)] = ii;
 
                 std::vector<int> bm(0);
                 bm.reserve(3*(submeshCell0[i + 1] - submeshCell0[i])*2);
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){
                     xM = 0;
                     yM = 0;
 //                    bm[1] = jj - submeshCell0[i];
                     for (int j = 0; j < 4; j++){
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int i0 = (int) std::round((xM - xMin[ii])/dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii])/dY[ii] - 0.5);
 
                     bm.push_back(i0 + j0*(2*ngc + nX[ii]));
                     bm.push_back(jj-submeshCell0[i]);
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
 
                 // enlarge TRIQUAD domain i
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){
 
                     double eps0 = std::min(dX[ii]/10, dY[ii]/10);
                     xM = 0;
                     yM = 0;
                     for (int j = 0; j < 4; j++){
                         nodes0[2*j] = xCoord[cells[jj][j]][0];
                         nodes0[2*j + 1] = xCoord[cells[jj][j]][1]-(yMax_temp-yMin_temp);
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1]-(yMax_temp-yMin_temp);
                     }
                     xM /= 4;
                     yM /= 4;
                     int cellFound = 0;
 
                     int l = submeshCell0[ii + nDomainsTRI];
                     bool foundCell = false;
                     while ((l < submeshCell0[ii + 1 + nDomainsTRI])&&((!foundCell)) ) {
                         int j = 0;
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++){
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
                         if ((std::abs(xM0 - xM) < eps0)&&(std::abs(yM0 - yM) < eps0)){
                             cells_add[i - (nDomainsQUAD + nDomainsTRI)][jj-submeshCell0[i]]= cells[l];
                             cellFound++;
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii])/dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii])/dY[ii] - 0.5);
 
                     bm.push_back(jj-submeshCell0[i] + submeshCell0[i + 1] - submeshCell0[i]);
                     bm.push_back(i0 + j0*(2*ngc + nX[ii]));
                     bm.push_back(-1); // note: -1 stands for interior update
 //                    assert(cellFound == 1 && "FastHybridMesh: Wrong number of cells found!");
                 }
                 bcMap_loc.push_back(bm);
                 break;
             }else if((std::fabs(yMin_temp - yMin[ii])<eps0)&&(std::fabs(xMax_temp-xMax[ii]+ngc*dX[ii])<eps0)&&
                      (std::fabs(xMin_temp-xMin[ii]-ngc*dX[ii])<eps0)){ //TRIQUAD domain i is below QUAD domain ii
                 // enlarge QUAD domain ii
 //                yMin[ii] = yMin_temp;
 //                nY[ii] += (yMax_temp - yMin_temp)/dY[ii];
                 assert((int)std::round((yMax_temp - yMin_temp)/dY[ii]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 TRIQUAD2QUAD[i - (nDomainsQUAD + nDomainsTRI)] = ii;
 
                 std::vector<int> bm(0);
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){
                     xM = 0;
                     yM = 0;
 //                    bm[1] = jj - submeshCell0[i];
                     for (int j = 0; j < 4; j++){
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int i0 = (int) std::round((xM - xMin[ii])/dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii])/dY[ii] - 0.5);
 
                     bm.push_back(i0 + j0*(2*ngc + nX[ii]));
                     bm.push_back(jj-submeshCell0[i]);
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
                 // enlarge TRIQUAD domain i
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){
 
                     double eps0 = std::min(dX[ii]/10, dY[ii]/10);
                     xM = 0;
                     yM = 0;
                     for (int j = 0; j < 4; j++){
                         nodes0[2*j] = xCoord[cells[jj][j]][0];
                         nodes0[2*j + 1] = xCoord[cells[jj][j]][1]+(yMax_temp-yMin_temp);
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1]+(yMax_temp-yMin_temp);
                     }
                     xM /= 4;
                     yM /= 4;
                     int cellFound = 0;
 
                     int l = submeshCell0[ii + nDomainsTRI];
                     bool foundCell = false;
                     while ((l < submeshCell0[ii + 1 + nDomainsTRI])&&((!foundCell)) ) {
                         int j = 0;
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++){
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
                         if ((std::abs(xM0 - xM) < eps0)&&(std::abs(yM0 - yM) < eps0)){
                             cells_add[i - (nDomainsQUAD + nDomainsTRI)][jj-submeshCell0[i]]= cells[l];
                             cellFound++;
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii])/dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii])/dY[ii] - 0.5);
 
                     bm.push_back(jj-submeshCell0[i] + submeshCell0[i + 1] - submeshCell0[i]);
                     bm.push_back(i0 + j0*(2*ngc + nX[ii]));
                     bm.push_back(-1); // note: -1 stands for interior update
 //                    assert(cellFound == 1 && "FastHybridMesh: Wrong number of cells found!");
                 }
                 bcMap_loc.push_back(bm);
                 break;
             }
         }
 
     }
 }
 
 void FastHybridMesh::addGhostCells2QUAD_2(std::vector<MeshType> &listPhysicalNameElements, std::vector<std::vector<std::vector<int>>> &cells_add,
                                         double & eps0, std::vector<std::vector<int>> &QUADQUAD2QUAD, std::vector<std::vector<int>> &bcMap_loc) {
     // for each QUADQUAD
 
     int i00 = nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT;
     for (int i = i00; i <  nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + nDomainsRT + nDomainsQUADQUAD;i++){
         QUADQUAD2QUAD[i - i00].resize(0);
         double xMin_temp = xCoord[cells[submeshCell0[i]][0]][0];
         double xMax_temp = xCoord[cells[submeshCell0[i]][0]][0];
         double yMin_temp = xCoord[cells[submeshCell0[i]][0]][1];
         double yMax_temp = xCoord[cells[submeshCell0[i]][0]][1];
 
         // find xMin, xMax, yMin, yMax for QUADQUAD
         for (int ii = submeshCell0[i]; ii < submeshCell0[i + 1]; ii++){
             for (int j = 0; j < 4; j++){
                 if (xMin_temp >  xCoord[cells[ii][j]][0]){
                     xMin_temp = xCoord[cells[ii][j]][0];
                 }
                 if (xMax_temp <  xCoord[cells[ii][j]][0]){
                     xMax_temp = xCoord[cells[ii][j]][0];
                 }
                 if (yMin_temp >  xCoord[cells[ii][j]][1]){
                     yMin_temp = xCoord[cells[ii][j]][1];
                 }
                 if (yMax_temp <  xCoord[cells[ii][j]][1]){
                     yMax_temp = xCoord[cells[ii][j]][1];
                 }
             }
         }
 //        int nCellsQUADQUAD = submeshCell0[i+1]-submeshCell0[i];
 //        cells_add[i - i00].resize(nCellsQUADQUAD, std::vector<int>(4));
         std::vector<double> nodes0(8);
 
         double xM, yM, xM0, yM0;
         for (int ii = 0; ii < nDomainsQUAD; ii++) {
             if ((std::fabs(xMax_temp-xMax[ii])<eps0)&&(std::fabs(yMax_temp-yMax[ii]+ngc*dY[ii])<eps0)&&
                 (std::fabs(yMin_temp-yMin[ii]-ngc*dY[ii])<eps0)){ //QUADQUAD domain i is right of QUAD domain ii
                 // enlarge QUAD domain ii and find Map from QUADQUAD to ghost cells of QUAD
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 QUADQUAD2QUAD[i - i00].push_back(ii);
 
                 // define QUADQUAD 2 QUAD map
                 std::vector<int> bm(0);
                 for (int j0 = 0; j0 < nY[ii]; j0++){
                     for (int i0 = 0; i0 < ngc; i0++) {
                         bm.push_back(i0 + nX[ii] + ngc + (j0 + ngc) * (2 * ngc + nX[ii]));
                         bm.push_back(i0 + ngc + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
 
 
                 // define QUAD 2 QUADQUAD map
                 for (int j0 = 0; j0 < nY[ii]; j0++){
                     for (int i0 = 0; i0 < ngc; i0++) {
                         bm.push_back(i0 + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD]));
                         bm.push_back(i0 + nX[ii] + (j0 + ngc) * (2 * ngc + nX[ii]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
                 bcMap_loc.push_back(bm);
             }
             if((std::fabs(xMin_temp - xMin[ii])<eps0)&&(std::fabs(yMax_temp-yMax[ii]+ngc*dY[ii])<eps0)&&
                      (std::fabs(yMin_temp-yMin[ii]-ngc*dY[ii])<eps0)){ //QUADQUAD domain i is left of QUAD domain ii
                 // enlarge QUAD domain ii
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 QUADQUAD2QUAD[i - i00].push_back(ii);
 
                 // define QUADQUAD 2 QUAD map
                 std::vector<int> bm(0);
                 for (int j0 = 0; j0 < nY[ii]; j0++){
                     for (int i0 = 0; i0 < ngc; i0++) {
                         bm.push_back(i0 + (j0 + ngc) * (2 * ngc + nX[ii]));
                         bm.push_back(i0 + nX[i - i00 + nDomainsQUAD] + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
 
 
                 // define QUAD 2 QUADQUAD map
                 for (int j0 = 0; j0 < nY[ii]; j0++){
                     for (int i0 = 0; i0 < ngc; i0++) {
                         bm.push_back(i0 + ngc + nX[i - i00 + nDomainsQUAD] + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD]));
                         bm.push_back(i0 + ngc + (j0 + ngc) * (2 * ngc + nX[ii]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
                 bcMap_loc.push_back(bm);
             }else if((std::fabs(yMax_temp - yMax[ii])<eps0)&&(std::fabs(xMax_temp-xMax[ii]+ngc*dX[ii])<eps0)&&
                      (std::fabs(xMin_temp-xMin[ii]-ngc*dX[ii])<eps0)){ //QUADQUAD domain i is on top of QUAD domain ii
                 // enlarge QUAD domain ii
                 assert((int) std::round((yMax_temp - yMin_temp)/dY[ii] )==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 QUADQUAD2QUAD[i - i00].push_back(ii);
 
                 // define QUADQUAD 2 QUAD map
                 std::vector<int> bm(0);
                 for (int i0 = 0; i0 < nX[ii]; i0++){
                     for (int j0 = 0; j0 < ngc; j0++) {
                         bm.push_back(i0 + ngc + (j0 + ngc + nY[ii]) * (2 * ngc + nX[ii]));
                         bm.push_back(i0 + ngc + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
 
 
                 // define QUAD 2 QUADQUAD map
                 for (int i0 = 0; i0 < nX[ii]; i0++){
                     for (int j0 = 0; j0 < ngc; j0++) {
                         bm.push_back(i0 + ngc + (j0) * (2 * ngc + nX[i - i00 + nDomainsQUAD]));
                         bm.push_back(i0 + ngc + (j0 + nY[ii]) * (2 * ngc + nX[ii]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
                 bcMap_loc.push_back(bm);
             }else if((std::fabs(yMin_temp - yMin[ii])<eps0)&&(std::fabs(xMax_temp-xMax[ii]+ngc*dX[ii])<eps0)&&
                      (std::fabs(xMin_temp-xMin[ii]-ngc*dX[ii])<eps0)){ //QUADQUAD domain i is below QUAD domain ii
                 // enlarge QUAD domain ii
 //                yMin[ii] = yMin_temp;
 //                nY[ii] += (yMax_temp - yMin_temp)/dY[ii];
                 assert((int)std::round((yMax_temp - yMin_temp)/dY[ii]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 QUADQUAD2QUAD[i - i00].push_back(ii);
 
                 // define QUADQUAD 2 QUAD map
                 std::vector<int> bm(0);
                 for (int i0 = 0; i0 < nX[ii]; i0++){
                     for (int j0 = 0; j0 < ngc; j0++) {
                         bm.push_back(i0 + ngc + (j0) * (2 * ngc + nX[ii]));
                         bm.push_back(i0 + ngc + (j0 + nY[i - i00 + nDomainsQUAD]) * (2 * ngc + nX[i - i00 + nDomainsQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
 
 
                 // define QUAD 2 QUADQUAD map
                 for (int i0 = 0; i0 < nX[ii]; i0++){
                     for (int j0 = 0; j0 < ngc; j0++) {
                         bm.push_back(i0 + ngc + (j0 + nY[i - i00 + nDomainsQUAD] + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD]));
                         bm.push_back(i0 + ngc + (j0 + ngc) * (2 * ngc + nX[ii]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
                 bcMap_loc.push_back(bm);
             }
         }
 
     }
 }
 
 void FastHybridMesh::addGhostCells2RT(std::vector<MeshType> &listPhysicalNameElements, std::vector<std::vector<std::vector<int>>> &cells_add,
                                         double & eps0, std::vector<std::vector<int>> &RT2QUADQUAD, std::vector<std::vector<int>> &bcMap_loc,
                                         std::vector<std::vector<int>> &QUADQUAD2QUAD, std::vector<std::vector<int>> &RT2QUAD) {
     // for each RT
     int i00 = nDomainsQUAD + nDomainsTRI + nDomainsTRIQUAD;
     for (int i = i00; i < i00 + nDomainsRT ;i++){
         double xMin_temp = xCoord[cells[submeshCell0[i]][0]][0];
         double xMax_temp = xCoord[cells[submeshCell0[i]][0]][0];
         double yMin_temp = xCoord[cells[submeshCell0[i]][0]][1];
         double yMax_temp = xCoord[cells[submeshCell0[i]][0]][1];
 
         // find xMin, xMax, yMin, yMax for RT
         for (int ii = submeshCell0[i]; ii < submeshCell0[i + 1]; ii++){
 
             for (int j = 0; j < 4; j++){
                 if (xMin_temp >  xCoord[cells[ii][j]][0]){
                     xMin_temp = xCoord[cells[ii][j]][0];
                 }
                 if (xMax_temp <  xCoord[cells[ii][j]][0]){
                     xMax_temp = xCoord[cells[ii][j]][0];
                 }
                 if (yMin_temp >  xCoord[cells[ii][j]][1]){
                     yMin_temp = xCoord[cells[ii][j]][1];
                 }
                 if (yMax_temp <  xCoord[cells[ii][j]][1]){
                     yMax_temp = xCoord[cells[ii][j]][1];
                 }
             }
         }
         int nCellsRT = submeshCell0[i+1]-submeshCell0[i];
         cells_add[i - i00 + nDomainsTRIQUAD].resize(0, std::vector<int>(4));
         std::vector<double> nodes0(8);
 
         double xM, yM, xM0, yM0;
         int cellFound = 0;
         // go through all QUADQUAD domains and check if they are direct neighbours
         for (int ii = 0; ii < nDomainsQUADQUAD; ii++) {
             if ((std::fabs(xMin_temp-xMin[ii + nDomainsQUAD])<eps0)&&(std::fabs(yMax_temp-yMax[ii + nDomainsQUAD]+ngc*dY[ii + nDomainsQUAD])<eps0)&&
                 (std::fabs(yMin_temp-yMin[ii + nDomainsQUAD]-ngc*dY[ii + nDomainsQUAD])<eps0)){ //QUADQUAD domain ii is right of RT domain i
                 // enlarge RT domain i and find Map from QUADQUAD to ghost cells of RT
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii + nDomainsQUAD]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 RT2QUADQUAD[i - i00].push_back(ii);
 
                 // RT 2 QUADQUAD map
                 std::vector<int> bm(0);
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){ // go through RT[i] cells
                     xM = 0;
                     yM = 0;
 
                     for (int j = 0; j < 4; j++){
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int i0 = (int) std::round((xM - xMin[ii + nDomainsQUAD])/dX[ii + nDomainsQUAD] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii + nDomainsQUAD])/dY[ii + nDomainsQUAD] - 0.5);
 
                     bm.push_back(i0 + j0*(2*ngc + nX[ii + nDomainsQUAD]));
                     bm.push_back(jj-submeshCell0[i]);
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
 
 
                 // enlarge RT domain i and find QUADQUAD 2 RT map
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){// go through RT[i] cells
 
                     eps0 = std::min(dX[ii + nDomainsQUAD]/10, dY[ii + nDomainsQUAD]/10);
 
                     xM = 0;
                     yM = 0;
                     for (int j = 0; j < 4; j++){
                         nodes0[2*j] = xCoord[cells[jj][j]][0] + ngc*dX[ii + nDomainsQUAD] ;
                         nodes0[2*j + 1] = xCoord[cells[jj][j]][1];
                         xM += xCoord[cells[jj][j]][0] + ngc*dX[ii + nDomainsQUAD];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
 
                     int l = submeshCell0[ii + i00 + nDomainsRT];
                     bool foundCell = false;
                     // go through all cells of QUADQUAD ii to find cell jj of RT
                     while ((l < submeshCell0[ii + i00 + nDomainsRT + 1])&&((!foundCell)) ) {
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++){
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
                         if ((std::abs(xM0 - xM) < eps0)&&(std::abs(yM0 - yM) < eps0)){
                             cells_add[i - i00 + nDomainsTRIQUAD].push_back(cells[l]);
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii + nDomainsQUAD])/dX[ii + nDomainsQUAD] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii + nDomainsQUAD])/dY[ii + nDomainsQUAD] - 0.5);
 
                     bm.push_back(submeshCell0[i + 1] - submeshCell0[i] + cells_add[i - i00 + nDomainsTRIQUAD].size()-1);
                     bm.push_back(i0 + j0*(2*ngc + nX[ii + nDomainsQUAD]));
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
                 bcMap_loc.push_back(bm);
             }else if((std::fabs(xMax_temp - xMax[ii + nDomainsQUAD])<eps0)&&(std::fabs(yMax_temp-yMax[ii + nDomainsQUAD]+ngc*dY[ii + nDomainsQUAD])<eps0)&&
                      (std::fabs(yMin_temp-yMin[ii + nDomainsQUAD]-ngc*dY[ii + nDomainsQUAD])<eps0)){ //QUADQUAD domain ii is left of RT domain i
                 // enlarge RT domain i and find Map from QUADQUAD to ghost cells of RT
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii + nDomainsQUAD]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 RT2QUADQUAD[i - i00].push_back(ii);
 
                 // RT 2 QUADQUAD map
                 std::vector<int> bm(0);
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){ // go through RT[i] cells
                     xM = 0;
                     yM = 0;
 
                     for (int j = 0; j < 4; j++){
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int i0 = (int) std::round((xM - xMin[ii + nDomainsQUAD])/dX[ii + nDomainsQUAD] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii + nDomainsQUAD])/dY[ii + nDomainsQUAD] - 0.5);
 
                     bm.push_back(i0 + j0*(2*ngc + nX[ii + nDomainsQUAD]));
                     bm.push_back(jj-submeshCell0[i]);
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
 
                 // enlarge RT domain i and find QUADQUAD 2 RT map
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){// go through RT[i] cells
 
                     eps0 = std::min(dX[ii + nDomainsQUAD]/10, dY[ii + nDomainsQUAD]/10);
 
                     xM = 0;
                     yM = 0;
                     for (int j = 0; j < 4; j++){
                         nodes0[2*j] = xCoord[cells[jj][j]][0] - ngc*dX[ii + nDomainsQUAD] ;
                         nodes0[2*j + 1] = xCoord[cells[jj][j]][1];
                         xM += xCoord[cells[jj][j]][0] - ngc*dX[ii + nDomainsQUAD];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
 
                     int l = submeshCell0[ii + i00 + nDomainsRT];
                     bool foundCell = false;
                     // go through all cells of QUADQUAD ii to find cell jj of RT
                     while ((l < submeshCell0[ii + i00 + nDomainsRT + 1])&&((!foundCell)) ) {
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++){
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
                         if ((std::abs(xM0 - xM) < eps0)&&(std::abs(yM0 - yM) < eps0)){
                             cells_add[i - i00 + nDomainsTRIQUAD].push_back(cells[l]);
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii + nDomainsQUAD])/dX[ii + nDomainsQUAD] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii + nDomainsQUAD])/dY[ii + nDomainsQUAD] - 0.5);
 
                     bm.push_back(submeshCell0[i + 1] - submeshCell0[i] + cells_add[i - i00 + nDomainsTRIQUAD].size()-1);
                     bm.push_back(i0 + j0*(2*ngc + nX[ii + nDomainsQUAD]));
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
                 bcMap_loc.push_back(bm);
             }else if((std::fabs(yMin_temp - yMin[ii + nDomainsQUAD])<eps0)&&(std::fabs(xMax_temp-xMax[ii + nDomainsQUAD]+ngc*dX[ii + nDomainsQUAD])<eps0)&&
                      (std::fabs(xMin_temp-xMin[ii + nDomainsQUAD]-ngc*dX[ii + nDomainsQUAD])<eps0)){ //QUADQUAD domain ii is on top of RT domain i
                 // enlarge RT domain i and find Map from QUADQUAD to ghost cells of RT
                 assert((int) std::round((yMax_temp - yMin_temp)/dY[ii + nDomainsQUAD] )==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 RT2QUADQUAD[i - i00].push_back(ii);
 
                 // RT 2 QUADQUAD map
                 std::vector<int> bm(0);
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){ // go through RT[i] cells
                     xM = 0;
                     yM = 0;
 
                     for (int j = 0; j < 4; j++){
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int i0 = (int) std::round((xM - xMin[ii + nDomainsQUAD])/dX[ii + nDomainsQUAD] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii + nDomainsQUAD])/dY[ii + nDomainsQUAD] - 0.5);
 
                     bm.push_back(i0 + j0*(2*ngc + nX[ii + nDomainsQUAD]));
                     bm.push_back(jj-submeshCell0[i]);
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
 
                 // enlarge RT domain i and find QUADQUAD 2 RT map
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){// go through RT[i] cells
 
                     eps0 = std::min(dX[ii + nDomainsQUAD]/10, dY[ii + nDomainsQUAD]/10);
 
                     xM = 0;
                     yM = 0;
                     for (int j = 0; j < 4; j++){
                         nodes0[2*j] = xCoord[cells[jj][j]][0] ;
                         nodes0[2*j + 1] = xCoord[cells[jj][j]][1] + ngc*dY[ii];
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1] + ngc*dY[ii];
                     }
                     xM /= 4;
                     yM /= 4;
 
                     int l = submeshCell0[ii + i00 + nDomainsRT];
                     bool foundCell = false;
                     // go through all cells of QUADQUAD ii to find cell jj of RT
                     while ((l < submeshCell0[ii + i00 + nDomainsRT + 1])&&((!foundCell)) ) {
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++){
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
                         if ((std::abs(xM0 - xM) < eps0)&&(std::abs(yM0 - yM) < eps0)){
                             cells_add[i - i00 + nDomainsTRIQUAD].push_back(cells[l]);
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii + nDomainsQUAD])/dX[ii + nDomainsQUAD] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii + nDomainsQUAD])/dY[ii + nDomainsQUAD] - 0.5);
 
                     bm.push_back(submeshCell0[i + 1] - submeshCell0[i] + cells_add[i - i00 + nDomainsTRIQUAD].size()-1);
                     bm.push_back(i0 + j0*(2*ngc + nX[ii + nDomainsQUAD]));
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
                 bcMap_loc.push_back(bm);
             }else if((std::fabs(yMax_temp - yMax[ii + nDomainsQUAD])<eps0)&&(std::fabs(xMax_temp-xMax[ii + nDomainsQUAD]+ngc*dX[ii + nDomainsQUAD])<eps0)&&
                      (std::fabs(xMin_temp-xMin[ii + nDomainsQUAD]-ngc*dX[ii + nDomainsQUAD])<eps0)){ //QUADQUAD domain ii is below of RT domain i
                 // enlarge RT domain i and find Map from QUADQUAD to ghost cells of RT
                 assert((int)std::round((yMax_temp - yMin_temp)/dY[ii + nDomainsQUAD]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 RT2QUADQUAD[i - i00].push_back(ii);
 
                 // RT 2 QUADQUAD map
                 std::vector<int> bm(0);
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){ // go through RT[i] cells
                     xM = 0;
                     yM = 0;
 
                     for (int j = 0; j < 4; j++){
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1];
                     }
                     xM /= 4;
                     yM /= 4;
                     int i0 = (int) std::round((xM - xMin[ii + nDomainsQUAD])/dX[ii + nDomainsQUAD] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii + nDomainsQUAD])/dY[ii + nDomainsQUAD] - 0.5);
 
                     bm.push_back(i0 + j0*(2*ngc + nX[ii + nDomainsQUAD]));
                     bm.push_back(jj-submeshCell0[i]);
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
 
                 // enlarge RT domain i and find QUADQUAD 2 RT map
                 for (int jj = submeshCell0[i]; jj < submeshCell0[i + 1]; jj++){// go through RT[i] cells
 
                     eps0 = std::min(dX[ii + nDomainsQUAD]/10, dY[ii + nDomainsQUAD]/10);
 
                     xM = 0;
                     yM = 0;
                     for (int j = 0; j < 4; j++){
                         nodes0[2*j] = xCoord[cells[jj][j]][0] ;
                         nodes0[2*j + 1] = xCoord[cells[jj][j]][1] - ngc*dY[ii + nDomainsQUAD];
                         xM += xCoord[cells[jj][j]][0];
                         yM += xCoord[cells[jj][j]][1] - ngc*dY[ii + nDomainsQUAD];
                     }
                     xM /= 4;
                     yM /= 4;
 
                     int l = submeshCell0[ii + i00 + nDomainsRT];
                     bool foundCell = false;
                     // go through all cells of QUADQUAD ii to find cell jj of RT
                     while ((l < submeshCell0[ii + i00 + nDomainsRT + 1])&&((!foundCell)) ) {
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++){
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
                         if ((std::abs(xM0 - xM) < eps0)&&(std::abs(yM0 - yM) < eps0)){
                             cells_add[i - i00 + nDomainsTRIQUAD].push_back(cells[l]);
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii + nDomainsQUAD])/dX[ii + nDomainsQUAD] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii + nDomainsQUAD])/dY[ii + nDomainsQUAD] - 0.5);
 
                     bm.push_back(submeshCell0[i + 1] - submeshCell0[i] + cells_add[i - i00 + nDomainsTRIQUAD].size()-1);
                     bm.push_back(i0 + j0*(2*ngc + nX[ii + nDomainsQUAD]));
                     bm.push_back(-1); // note: -1 stands for interior update
                 }
                 bcMap_loc.push_back(bm);
             }
         }
         // find QUAD that is between the QUADQUADs
         if (RT2QUADQUAD[i-i00].size() == 2){
             int iQUAD;
             if (QUADQUAD2QUAD[RT2QUADQUAD[i-i00][0]][0] == QUADQUAD2QUAD[RT2QUADQUAD[i-i00][1]][0]){
                 iQUAD = QUADQUAD2QUAD[RT2QUADQUAD[i-i00][0]][0];
             }else if (QUADQUAD2QUAD[RT2QUADQUAD[i-i00][0]][0] == QUADQUAD2QUAD[RT2QUADQUAD[i-i00][1]][1]){
                 iQUAD = QUADQUAD2QUAD[RT2QUADQUAD[i-i00][0]][0];
             }else{
                 iQUAD = QUADQUAD2QUAD[RT2QUADQUAD[i-i00][0]][1];
             }
             RT2QUAD[i-i00].push_back(iQUAD);
 
             int imin, imax, jmin, jmax;
             // find out how QUAD and RT are positioned
             if (std::fabs(xMin_temp - (xMax[iQUAD]-ngc*dX[iQUAD])) < eps0&&std::fabs(yMin_temp - (yMax[iQUAD]-ngc*dX[iQUAD]))< eps0){ // QUAD is south west of RT
                 imin = nX[iQUAD];
                 imax = nX[iQUAD] + ngc;
                 jmin = nY[iQUAD];
                 jmax = nY[iQUAD] + ngc;
             }else if (std::fabs(xMax_temp - (xMin[iQUAD]+ngc*dX[iQUAD]))< eps0&&std::fabs(yMin_temp - (yMax[iQUAD]-ngc*dX[iQUAD]))< eps0){ // QUAD is south east of RT
                 imin = ngc;
                 imax = 2*ngc;
                 jmin = nY[iQUAD];
                 jmax = nY[iQUAD] + ngc;
             }else if (std::fabs(xMax_temp - (xMin[iQUAD]+ngc*dX[iQUAD]))< eps0&&std::fabs(yMax_temp - (yMin[iQUAD]+ngc*dX[iQUAD]))< eps0){ // QUAD is north east of RT
                 imin = ngc;
                 imax = 2*ngc;
                 jmin = ngc;
                 jmax = 2*ngc;
             }else{// QUAD is north west of RT
                 imin = nX[iQUAD];
                 imax = nX[iQUAD] + ngc;
                 jmin = ngc;
                 jmax = 2*ngc;
             }
             std::vector<int> bm(0);
             for (int ii = imin; ii < imax; ii++){
                 for (int jj = jmin; jj < jmax; jj++) {
 
                     double eps0 = std::min(dX[iQUAD] / 10, dY[iQUAD] / 10);
                     xM = (ii + 0.5) * dX[iQUAD] + xMin[iQUAD];
                     yM = (jj + 0.5) * dY[iQUAD] + yMin[iQUAD];
                     int cellFound = 0;
 
                     int l = submeshCell0[iQUAD + nDomainsTRI];
                     bool foundCell = false;
                     while ((l < submeshCell0[iQUAD + 1 + nDomainsTRI]) && ((!foundCell))) {
                         int j = 0;
                         xM0 = 0;
                         yM0 = 0;
                         for (int j = 0; j < 4; j++) {
                             xM0 += xCoord[cells[l][j]][0];
                             yM0 += xCoord[cells[l][j]][1];
                         }
 
                         xM0 /= 4;
                         yM0 /= 4;
 
 //                        std::cout << std::abs(xM0 - xM) << ", " << std::abs(yM0 - yM) << "\n";
                         if ((std::abs(xM0 - xM) < eps0) && (std::abs(yM0 - yM) < eps0)) {
                             cells_add[i - i00 + nDomainsTRIQUAD].push_back(cells[l]);
                             cellFound++;
                             foundCell = true;
                         }
                         l++;
                     }
                     int i0 = (int) std::round((xM - xMin[ii]) / dX[ii] - 0.5);
                     int j0 = (int) std::round((yM - yMin[ii]) / dY[ii] - 0.5);
 
                     bm.push_back(submeshCell0[i + 1] - submeshCell0[i] + cells_add[i - i00 + nDomainsTRIQUAD].size()-1);
                     bm.push_back(ii + jj*(2*ngc + nX[iQUAD]));
                     bm.push_back(-1);
                 }
             }
             bcMap_loc.push_back(bm);
         }
     }
 }
 
 
 void FastHybridMesh::addGhostCells2RQ(std::vector<MeshType> &listPhysicalNameElements, std::vector<std::vector<std::vector<int>>> &cells_add,
 double & eps0, std::vector<std::vector<int>> &RQ2QUAD, std::vector<std::vector<int>> &bcMap_loc){
 // for each RQ
     int i00 = nDomainsQUAD + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsQUADQUAD;
     for (int i = i00; i < i00 + nDomainsRQ ;i++){
         double xMin_temp = xCoord[cells[submeshCell0[i]][0]][0];
         double xMax_temp = xCoord[cells[submeshCell0[i]][0]][0];
         double yMin_temp = xCoord[cells[submeshCell0[i]][0]][1];
         double yMax_temp = xCoord[cells[submeshCell0[i]][0]][1];
 
         // find xMin, xMax, yMin, yMax for RQ
         for (int ii = submeshCell0[i]; ii < submeshCell0[i + 1]; ii++){
 
             for (int j = 0; j < 4; j++){
                 if (xMin_temp >  xCoord[cells[ii][j]][0]){
                     xMin_temp = xCoord[cells[ii][j]][0];
                 }
                 if (xMax_temp <  xCoord[cells[ii][j]][0]){
                     xMax_temp = xCoord[cells[ii][j]][0];
                 }
                 if (yMin_temp >  xCoord[cells[ii][j]][1]){
                     yMin_temp = xCoord[cells[ii][j]][1];
                 }
                 if (yMax_temp <  xCoord[cells[ii][j]][1]){
                     yMax_temp = xCoord[cells[ii][j]][1];
                 }
             }
         }
         int nCellsRQ = submeshCell0[i+1]-submeshCell0[i];
 //        cells_add[i - i00 + nDomainsTRIQUAD + nDomainsRT].resize(nCellsRQ, std::vector<int>(4));
         std::vector<double> nodes0(8);
 
         double xM, yM, xM0, yM0;
         // go through all QUADQUAD domains and check if they are direct neighbours
         for (int ii = 0; ii < nDomainsQUADQUAD; ii++) {
             if ((std::fabs(xMin_temp-xMin[ii + nDomainsQUAD])<eps0)&&(std::fabs(yMax_temp-yMax[ii + nDomainsQUAD]+ngc*dY[ii + nDomainsQUAD])<eps0)&&
                 (std::fabs(yMin_temp-yMin[ii + nDomainsQUAD]-ngc*dY[ii + nDomainsQUAD])<eps0)){ //QUADQUAD domain ii is right of RQ domain i
                 // find Map from QUADQUAD to ghost cells of RQ
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii + nDomainsQUAD]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 RQ2QUAD[i - i00].push_back(ii);
 
                 // RQ 2 QUADQUAD map
                 std::vector<int> bm(0);
                 for (int j0 = 0; j0 < nY[ii + nDomainsQUAD]; j0++){
                     for (int i0 = 0; i0 < ngc; i0++) {
                         bm.push_back(i0 + (j0 + ngc) * (2 * ngc + nX[ii + nDomainsQUAD]));
                         bm.push_back(i0 + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD] + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
 
                 // define QUADQUAD 2 RQ map
                 for (int j0 = 0; j0 < nY[ii + nDomainsQUAD]; j0++){
                     for (int i0 = 0; i0 < ngc; i0++) {
                         bm.push_back(i0 + ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD] + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]));
                         bm.push_back(i0 + ngc + (j0 + ngc) * (2 * ngc + nX[ii + nDomainsQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
                 bcMap_loc.push_back(bm);
             }else if((std::fabs(xMax_temp - xMax[ii + nDomainsQUAD])<eps0)&&(std::fabs(yMax_temp-yMax[ii + nDomainsQUAD]+ngc*dY[ii + nDomainsQUAD])<eps0)&&
                      (std::fabs(yMin_temp-yMin[ii + nDomainsQUAD]-ngc*dY[ii + nDomainsQUAD])<eps0)){ //QUADQUAD domain ii is left of RQ domain i
                 // find Map from QUADQUAD to ghost cells of RQ
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii + nDomainsQUAD]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 RQ2QUAD[i - i00].push_back(ii);
 
                 // RQ 2 QUADQUAD map
                 std::vector<int> bm(0);
                 for (int j0 = 0; j0 < nY[ii + nDomainsQUAD]; j0++){
                     for (int i0 = 0; i0 < ngc; i0++) {
                         bm.push_back(i0 + ngc + nX[ii + nDomainsQUAD] + (j0 + ngc) * (2 * ngc + nX[ii + nDomainsQUAD]));
                         bm.push_back(i0 + ngc + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
 
                 // define QUADQUAD 2 RQ map
                 for (int j0 = 0; j0 < nY[ii + nDomainsQUAD]; j0++){
                     for (int i0 = 0; i0 < ngc; i0++) {
                         bm.push_back(i0 + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]));
                         bm.push_back(i0 + nX[ii + nDomainsQUAD] + (j0 + ngc) * (2 * ngc + nX[ii + nDomainsQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
                 bcMap_loc.push_back(bm);
             }else if((std::fabs(yMin_temp - yMin[ii + nDomainsQUAD])<eps0)&&(std::fabs(xMax_temp-xMax[ii + nDomainsQUAD]+ngc*dX[ii + nDomainsQUAD])<eps0)&&
                      (std::fabs(xMin_temp-xMin[ii + nDomainsQUAD]-ngc*dX[ii + nDomainsQUAD])<eps0)){ //QUADQUAD domain ii is on top of RQ domain i
                 // find Map from QUADQUAD to ghost cells of RQ
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii + nDomainsQUAD]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 RQ2QUAD[i - i00].push_back(ii);
 
                 // RQ 2 QUADQUAD map
                 std::vector<int> bm(0);
                 for (int i0 = 0; i0 < nX[ii + nDomainsQUAD]; i0++){
                     for (int j0 = 0; j0 < ngc; j0++) {
                         bm.push_back(i0 + ngc + (j0) * (2 * ngc + nX[ii + nDomainsQUAD]));
                         bm.push_back(i0 + ngc + (j0 + nY[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]) * (2 * ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
 
                 // define QUADQUAD 2 RQ map
                 for (int i0 = 0; i0 < nX[ii + nDomainsQUAD]; i0++){
                     for (int j0 = 0; j0 < ngc; j0++) {
                         bm.push_back(i0 + ngc + (j0 + ngc + nY[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]) * (2 * ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]));
                         bm.push_back(i0 + ngc + (j0 + ngc) * (2 * ngc + nX[ii + nDomainsQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
                 bcMap_loc.push_back(bm);
             }else if((std::fabs(yMax_temp - yMax[ii + nDomainsQUAD])<eps0)&&(std::fabs(xMax_temp-xMax[ii + nDomainsQUAD]+ngc*dX[ii + nDomainsQUAD])<eps0)&&
                      (std::fabs(xMin_temp-xMin[ii + nDomainsQUAD]-ngc*dX[ii + nDomainsQUAD])<eps0)){ //QUADQUAD domain ii is below of RQ domain i
                 // find Map from QUADQUAD to ghost cells of RQ
                 assert((int) std::round((xMax_temp - xMin_temp)/dX[ii + nDomainsQUAD]) ==ngc && "FastHybridMesh: Wrong number of ghost cells in mesh!");
                 RQ2QUAD[i - i00].push_back(ii);
 
                 // RQ 2 QUADQUAD map
                 std::vector<int> bm(0);
                 for (int i0 = 0; i0 < nX[ii + nDomainsQUAD]; i0++){
                     for (int j0 = 0; j0 < ngc; j0++) {
                         bm.push_back(i0 + ngc + (j0 + ngc + nY[ii + nDomainsQUAD]) * (2 * ngc + nX[ii + nDomainsQUAD]));
                         bm.push_back(i0 + ngc + (j0 + ngc) * (2 * ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
 
                 // define QUADQUAD 2 RQ map
                 for (int i0 = 0; i0 < nX[ii + nDomainsQUAD]; i0++){
                     for (int j0 = 0; j0 < ngc; j0++) {
                         bm.push_back(i0 + ngc + (j0) * (2 * ngc + nX[i - i00 + nDomainsQUAD + nDomainsQUADQUAD]));
                         bm.push_back(i0 + ngc + (j0 + nY[ii + nDomainsQUAD]) * (2 * ngc + nX[ii + nDomainsQUAD]));
                         bm.push_back(-1); // note: -1 stands for interior update
                     }
                 }
                 bcMap_loc.push_back(bm);
             }
         }
 
     }
 }
 
 void FastHybridMesh::splitCellsForRBF(std::vector<std::vector<std::vector<int>>> &cells_add){
 
     int tmpCell;
     int n0, n1;
 
     // Split cells of RT (in cells)
     gmsh::logger::write("Split cells of RT (in cells).");
     for (int i = 0; i < nDomainsRT; i++){
         n0 = submeshCell0[i + nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD];
         n1 = submeshCell0[i + nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + 1];
         for (int j = n0; j < n1; j++){
             tmpCell = cells[j][3];
             cells[j].pop_back();
             cells.insert(cells.begin() + submeshCell0[i + nDomainsTRI + nDomainsQUAD + nDomainsTRIQUAD + 1], std::vector<int>{cells[j][2], tmpCell, cells[j][0]});
             for (int l = i + nDomainsQUAD + nDomainsTRI + nDomainsTRIQUAD + 1; l < submeshCell0.size(); l++){
                 submeshCell0[l] += 1;
             }
         }
     }
 
     // Split cells of TRIQUAD (in cells)
     gmsh::logger::write("Split cells of TRIQUAD (in cells).");
     for (int i = 0; i < nDomainsTRIQUAD; i++){
         n0 = submeshCell0[i + nDomainsTRI + nDomainsQUAD];
         n1 = submeshCell0[i + nDomainsTRI + nDomainsQUAD + 1];
         for (int j = n0; j < n1; j++){
             tmpCell = cells[j][3];
             cells[j].pop_back();
             cells.insert(cells.begin() + submeshCell0[i + nDomainsTRI + nDomainsQUAD + 1], std::vector<int>{cells[j][2], tmpCell, cells[j][0]});
             for (int l = i + nDomainsQUAD + nDomainsTRI + 1; l < submeshCell0.size(); l++){
                 submeshCell0[l] += 1;//submeshCell0[i + nDomainsTRI + nDomainsQUAD + 1] - submeshCell0[i + nDomainsTRI + nDomainsQUAD];
             }
         }
     }
 
     // Split cells_add cells (for TRIQUAD) and create map of splitting
     gmsh::logger::write("Split cells_add cells (for TRIQUAD) and create map of splitting.");
     int nCellsAdd;
 
     for (int i = nDomainsQUAD + nDomainsTRI; i <  nDomainsQUAD + nDomainsTRI + nDomainsTRIQUAD;i++){
         nCellsAdd = cells_add[i - (nDomainsQUAD + nDomainsTRI)].size();
         for (int j = 0;j < nCellsAdd; j++){
             tmpCell = cells_add[i - (nDomainsQUAD + nDomainsTRI)][j][3];
             cells_add[i - (nDomainsQUAD + nDomainsTRI)][j].pop_back();
             cells_add[i - (nDomainsQUAD + nDomainsTRI)].push_back(std::vector<int>{cells_add[i - (nDomainsQUAD + nDomainsTRI)][j][2], tmpCell, cells_add[i - (nDomainsQUAD + nDomainsTRI)][j][0]});
         }
     }
 
 
 
     // Split cells_add cells (for RT) and create map of splitting
     gmsh::logger::write("Split cells_add cells (for RT) and create map of splitting.");
 
     for (int i = nDomainsQUAD + nDomainsTRI + nDomainsTRIQUAD; i <  nDomainsQUAD + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT;i++){
         nCellsAdd = cells_add[i - (nDomainsQUAD + nDomainsTRI)].size();
         for (int j = 0;j < nCellsAdd; j++){
             tmpCell = cells_add[i - (nDomainsQUAD + nDomainsTRI)][j][3];
             cells_add[i - (nDomainsQUAD + nDomainsTRI)][j].pop_back();
             cells_add[i - (nDomainsQUAD + nDomainsTRI)].push_back(std::vector<int>{cells_add[i
                                                                                              - (nDomainsQUAD + nDomainsTRI)][j][2], tmpCell, cells_add[i - (nDomainsQUAD + nDomainsTRI)][j][0]});
         }
     }
 }
 
 void FastHybridMesh::setGlobalBcMap(std::vector<std::vector<int>> &bcMap_loc, std::vector<int> &TRIQUAD2QUAD, std::vector<std::vector<int>> &QUADQUAD2QUAD,
         std::vector<std::vector<int>> &RQ2QUADQUAD, std::vector<std::vector<int>> &RT2QUADQUAD, std::vector<std::vector<int>> &RT2QUAD){
     int i0, i1;
     // TRIQUAD and QUAD
     for (int i = 0; i < nDomainsTRIQUAD; i++){
         int n = bcMap_loc[i].size()/6;
 
         //TRIQUAD 2 QUAD
         i0 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsRQ + nDomainsQUADQUAD;
         i1 = nDomainsTRI;
         for (int j = 0; j < n; j++){
             bcMap2.push_back(bcMap_loc[i][3*j + 0] + submeshCell0[i0 + TRIQUAD2QUAD[i]]);
             bcMap2.push_back(bcMap_loc[i][3*j + 1] + submeshCell0[i1 + i]);
             bcMap2.push_back(-1);
         }
 
 
         // QUAD 2 TRIQUAD
         for (int j = 0; j < n; j++){
             bcMap.push_back(bcMap_loc[i][3*(j+n) + 0] + submeshCell0[i1 + i] + n);
             bcMap.push_back(bcMap_loc[i][3*(j+n) + 1] + submeshCell0[i0 + TRIQUAD2QUAD[i]]);
             bcMap.push_back(-1);
             bcMap.push_back(bcMap_loc[i][3*(j+n) + 0] + submeshCell0[i1 + i] + 2*n);
             bcMap.push_back(bcMap_loc[i][3*(j+n) + 1] + submeshCell0[i0 + TRIQUAD2QUAD[i]]);
             bcMap.push_back(-1);
         }
     }
 
 //    for (int i = 0; i < nDomainsTRIQUAD; i++){
 //        int n = bcMap_loc[i].size()/6;
 //
 //        //TRIQUAD 2 QUAD
 //        for (int j = 0; j < n; j++){
 //            bcMap2.push_back(bcMap_loc[i][3*j + 0] + submeshCell0[i0 + TRIQUAD2QUAD[i]]);
 //            bcMap2.push_back(bcMap_loc[i][3*j + 1] + submeshCell0[i1 + i] + n);
 //            bcMap2.push_back(-1);
 //        }
 //    }
 
     // QUADQUAD and QUAD
     for (int i = 0; i < nDomainsQUADQUAD; i++){
         int n = bcMap_loc[2*i + nDomainsTRIQUAD].size()/6;
 
         //QUADQUAD 2 QUAD
         i0 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsRQ + nDomainsQUADQUAD;
         i1 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT;
         for (int j = 0; j < n; j++){
             bcMap.push_back(bcMap_loc[2*i + nDomainsTRIQUAD][3*j + 0] + submeshCell0[i0 + QUADQUAD2QUAD[i][0]]);
             bcMap.push_back(bcMap_loc[2*i + nDomainsTRIQUAD][3*j + 1] + submeshCell0[i1 + i]);
             bcMap.push_back(-1);
         }
 
         // QUAD 2 QUADQUAD
         for (int j = 0; j < n; j++){
             bcMap.push_back(bcMap_loc[2*i + nDomainsTRIQUAD][3*(j+n) + 0] + submeshCell0[i1 + i]);
             bcMap.push_back(bcMap_loc[2*i + nDomainsTRIQUAD][3*(j+n) + 1] + submeshCell0[i0 + QUADQUAD2QUAD[i][0]]);
             bcMap.push_back(-1);
         }
         n = bcMap_loc[2*i + 1 + nDomainsTRIQUAD].size()/6;
 
         //QUADQUAD 2 QUAD
         for (int j = 0; j < n; j++){
             bcMap.push_back(bcMap_loc[2*i + 1 + nDomainsTRIQUAD][3*j + 0] + submeshCell0[i0 + QUADQUAD2QUAD[i][1]]);
             bcMap.push_back(bcMap_loc[2*i + 1 + nDomainsTRIQUAD][3*j + 1] + submeshCell0[i1 + i]);
             bcMap.push_back(-1);
         }
 
         // QUAD 2 QUADQUAD
         for (int j = 0; j < n; j++){
             bcMap.push_back(bcMap_loc[2*i + 1 + nDomainsTRIQUAD][3*(j+n) + 0] + submeshCell0[i1 + i]);
             bcMap.push_back(bcMap_loc[2*i + 1 + nDomainsTRIQUAD][3*(j+n) + 1] + submeshCell0[i0 + QUADQUAD2QUAD[i][1]]);
             bcMap.push_back(-1);
         }
     }
 
     // RQ and QUADQUAD
     int i0bcMapRQ = bcMap_loc.size() - nDomainsRQ*4;
     for (int i = 0; i < nDomainsRQ; i++){
         for (int jQQ = 0; jQQ < 4; jQQ++) {
             int n = bcMap_loc[4 * i + jQQ + i0bcMapRQ].size() / 6;
 
             //RQ 2 QUADQUAD
             i0 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT;
             i1 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsQUADQUAD;
             for (int j = 0; j < n; j++) {
                 bcMap.push_back(bcMap_loc[4 * i + jQQ + i0bcMapRQ][3 * j + 0] + submeshCell0[i0 + RQ2QUADQUAD[i][jQQ]]);
                 bcMap.push_back(bcMap_loc[4 * i + jQQ + i0bcMapRQ][3 * j + 1] + submeshCell0[i1 + i]);
                 bcMap.push_back(-1);
             }
 
             // QUADQUAD 2 RQ
             for (int j = 0; j < n; j++) {
                 bcMap.push_back(bcMap_loc[4 * i + jQQ + i0bcMapRQ][3 * (j + n) + 0] + submeshCell0[i1 + i]);
                 bcMap.push_back(bcMap_loc[4 * i + jQQ + i0bcMapRQ][3 * (j + n) + 1] + submeshCell0[i0 + RQ2QUADQUAD[i][jQQ]]);
                 bcMap.push_back(-1);
             }
         }
     }
 
     // RT and QUADQUAD
     int indRT = -1;
     for (int i = 0; i < nDomainsRT; i++){
         int n0 = 0;
         for (int l = 0; l < RT2QUADQUAD[i].size(); l++){
             n0 += bcMap_loc[indRT + 1 + l + nDomainsTRIQUAD + 2*nDomainsQUADQUAD].size()/6;
         }
         for (int l = 0; l < RT2QUAD[i].size(); l++){
             n0 += bcMap_loc[indRT + RT2QUADQUAD[i].size() + l + nDomainsTRIQUAD + 2*nDomainsQUADQUAD].size()/6;
         }
 
         for (int l = 0; l < RT2QUADQUAD[i].size(); l++){
             indRT++;
             int n = bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD].size()/6;
 
             //RT 2 QUADQUAD
             i0 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT;
             i1 = nDomainsTRI + nDomainsTRIQUAD;
             for (int j = 0; j < n; j++){
                 bcMap2.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*j + 0] + submeshCell0[i0 + RT2QUADQUAD[i][l]]);
                 bcMap2.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*j + 1] + submeshCell0[i1 + i]);
                 bcMap2.push_back(-1);
             }
 
 
             // QUAD 2 RT
             for (int j = 0; j < n; j++){
                 bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*(j+n) + 0] + submeshCell0[i1 + i] + n);
                 bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*(j+n) + 1] + submeshCell0[i0 + RT2QUADQUAD[i][l]]);
                 bcMap.push_back(-1);
                 bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*(j+n) + 0] + submeshCell0[i1 + i] + n + n0);
                 bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*(j+n) + 1] + submeshCell0[i0 + RT2QUADQUAD[i][l]]);
                 bcMap.push_back(-1);
             }
         }
 
         // QUAD 2 RT
         for (int l = 0; l < RT2QUAD[i].size(); l++){
             indRT++;
             int n = bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD].size()/3;
             int n1 = 0;
             for (int l = 0; l < RT2QUADQUAD[i].size(); l++){
                 n1 += bcMap_loc[indRT + 1 + l + nDomainsTRIQUAD + 2*nDomainsQUADQUAD].size()/6;
             }
             //RT 2 QUADQUAD
             i0 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsQUADQUAD;
             i1 = nDomainsTRI + nDomainsTRIQUAD;
             for (int j = 0; j < n; j++){
                 bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*(j) + 0] + submeshCell0[i1 + i] + n);
                 bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*(j) + 1] + submeshCell0[i0 + RT2QUAD[i][l]]);
                 bcMap.push_back(-1);
                 bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*(j) + 0] + submeshCell0[i1 + i] + n + n0);
                 bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*(j) + 1] + submeshCell0[i0 + RT2QUAD[i][l]]);
                 bcMap.push_back(-1);
             }
         }
     }
 
     for (int i = 0; i < nDomainsTRIQUAD; i++){
         int n = bcMap_loc[i].size()/6;
         i0 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + nDomainsRQ + nDomainsQUADQUAD;
         i1 = nDomainsTRI;
         //TRIQUAD 2 QUAD
         for (int j = 0; j < n; j++){
             bcMap2.push_back(bcMap_loc[i][3*j + 0] + submeshCell0[i0 + TRIQUAD2QUAD[i]]);
             bcMap2.push_back(bcMap_loc[i][3*j + 1] + submeshCell0[i1 + i] + n);
             bcMap2.push_back(-1);
         }
     }
 
     indRT = -1;
     for (int i = 0; i < nDomainsRT; i++){
         int n0 = 0;
         for (int l = 0; l < RT2QUADQUAD[i].size(); l++){
             n0 += bcMap_loc[indRT + 1 + l + nDomainsTRIQUAD + 2*nDomainsQUADQUAD].size()/6;
         }
         for (int l = 0; l < RT2QUAD[i].size(); l++){
             n0 += bcMap_loc[indRT + RT2QUADQUAD[i].size() + l + nDomainsTRIQUAD + 2*nDomainsQUADQUAD].size()/6;
         }
         for (int l = 0; l < RT2QUADQUAD[i].size(); l++){
             indRT++;
             int n = bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD].size()/6;
 
             //RT 2 QUADQUAD
             i0 = nDomainsTRI + nDomainsTRIQUAD + nDomainsRT;
             i1 = nDomainsTRI + nDomainsTRIQUAD;
             for (int j = 0; j < n; j++){
                 bcMap2.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*j + 0] + submeshCell0[i0 + RT2QUADQUAD[i][l]]);
                 bcMap2.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2*nDomainsQUADQUAD][3*j + 1] + submeshCell0[i1 + i] + n);
                 bcMap2.push_back(-1);
             }
         }
     }
 }
 
 std::vector<std::vector<double>>  FastHybridMesh::getXCoordOfCell(int i){
-    std::vector<std::vector<double>> ret;
-    int n = submeshCell0.size();
-    int j = 1;
-    while (submeshCell0[j] < i ){
-        j++;
-    }
-
-    if (i < nDomainsTRI + nDomainsTRIQUAD + nDomainsRT){
-        ret.resize(3,std::vector<double>(2,0));
-        std::vector<int> cell_i = cells[i];
-        ret[0][0] = xCoord[cell_i[0]][0];
-        ret[0][1] = xCoord[cell_i[0]][1];
-        ret[1][0] = xCoord[cell_i[1]][0];
-        ret[1][1] = xCoord[cell_i[1]][1];
-        ret[2][0] = xCoord[cell_i[2]][0];
-        ret[2][1] = xCoord[cell_i[2]][1];
-    }
-    else{
-        j = j - nDomainsTRI + nDomainsTRIQUAD + nDomainsRT;
-        int j0 = floor(((double)(i))/(nX[j] + 2 * ngc));
-        int i0 = j - j0*(nX[j] + 2 * ngc);
-    }
+//    std::vector<std::vector<double>> ret;
+//    int n = submeshCell0.size();
+//    int j = 1;
+//    while (submeshCell0[j] < i ){
+//        j++;
+//    }
+//
+//    if (j < nDomainsTRI + nDomainsTRIQUAD + nDomainsRT){
+//        ret.resize(3,std::vector<double>(2,0));
+//        std::vector<int> cell_i = cells[i];
+//        ret[0][0] = xCoord[cell_i[0]][0];
+//        ret[0][1] = xCoord[cell_i[0]][1];
+//        ret[1][0] = xCoord[cell_i[1]][0];
+//        ret[1][1] = xCoord[cell_i[1]][1];
+//        ret[2][0] = xCoord[cell_i[2]][0];
+//        ret[2][1] = xCoord[cell_i[2]][1];
+//    }else{
+//        j = j - nDomainsTRI + nDomainsTRIQUAD + nDomainsRT;
+//        int j0 = floor(((double)(i))/(nX[j] + 2 * ngc));
+//        int i0 = j - j0*(nX[j] + 2 * ngc);
+//    }
+    std::vector<std::vector<double>> ret(3,std::vector<double>(2,0));
+    std::vector<int> cell_i = cells[i];
+    ret[0][0] = xCoord[cell_i[0]][0];
+    ret[0][1] = xCoord[cell_i[0]][1];
+    ret[1][0] = xCoord[cell_i[1]][0];
+    ret[1][1] = xCoord[cell_i[1]][1];
+    ret[2][0] = xCoord[cell_i[2]][0];
+    ret[2][1] = xCoord[cell_i[2]][1];
+    return ret;
 }
 
 # undef NODE_NUM
\ No newline at end of file