FastHybridMesh.cpp
031392caa76amaster

FastHybridMesh.cpp
View Options

			//
			// 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] + 2ngc)(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] + 2ngc)(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]+2ngc)(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 = ngcngc2;
			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 = ngcngc2;
			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] + 2ngc);
			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] + 2ngc);
			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] + 2ngc);
			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 = ngcngc2;
			// 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] + 2ngc));
			// }
			// }
			// }
			//
			// 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] + 2ngc));
			// }
			// }
			// }
			//
			// 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] + 2ngc));
			// }
			// }
			// }

			// 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] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + iii];
			// bm[1] = i + ngc + (2ngc - 1 - jcell)(nX[iii] + 2*ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + iii];
			// bm[2] = edgesAtBC.size()-1;
			bcMap.push_back(i + ngc + (jcell)(nX[iii] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
			if (edgeInfo[1] == ABSORBINGBC){
			bcMap.push_back(i + ngc + (ngc)(nX[iii] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
			}else{
			bcMap.push_back(i + ngc + (2ngc - 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] + 2ngcdY[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] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
			if (edgeInfo[1] == ABSORBINGBC){
			bcMap.push_back(i + ngc + ( - 1 + nY[iii] + ngc)(nX[iii] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
			} else{
			bcMap.push_back(i + ngc + (-jcell - 1 + nY[iii] + ngc)(nX[iii] + 2ngc) + 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] + 2ngcdY[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] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
			if (edgeInfo[1] == ABSORBINGBC){
			bcMap.push_back(ngc + (j + ngc)(nX[iii] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
			}else{
			bcMap.push_back(2ngc - 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] + 2ngcdX[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] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
			if (edgeInfo[1] == ABSORBINGBC){
			bcMap.push_back(nX[iii] + ngc - 1 + (j + ngc)(nX[iii] + 2ngc) + submeshCell0[nDomainsTRI + nDomainsTRIQUAD + nDomainsRT + iii]);
			}else{
			bcMap.push_back(nX[iii] + ngc - icell - 1 + (j + ngc)(nX[iii] + 2ngc) + 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] + 2ngcdX[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(2nn);
			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]- (iidX[i] + xMin[i]+ngcdX[i]))< eps0) &&
			// (std::abs(xCoord[k][1]- (jjdY[i] + yMin[i]+ngcdY[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]- (iidX[i + nDomainsQUAD] + xMin[i + nDomainsQUAD]+ngcdX[i + nDomainsQUAD]))< eps0) &&
			// (std::abs(xCoord[k][1]- (jjdY[i + nDomainsQUAD] + yMin[i + nDomainsQUAD]+ngcdY[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]- (iidX[i + nDomainsQUAD + nDomainsQUADQUAD] + xMin[i + nDomainsQUAD + nDomainsQUADQUAD]+ngcdX[i + nDomainsQUAD + nDomainsQUADQUAD]))< eps0) &&
			// (std::abs(xCoord[k][1]- (jjdY[i + nDomainsQUAD + nDomainsQUADQUAD] + yMin[i + nDomainsQUAD + nDomainsQUADQUAD]+ngcdY[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]- (iidX[i + nDomainsQUAD] + xMin[i + nDomainsQUAD]+ngcdX[i + nDomainsQUAD]))< eps0) &&
			(std::abs(xCoord[k][1]- (jjdY[i + nDomainsQUAD] + yMin[i + nDomainsQUAD]+ngcdY[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]- (iidX[i + nDomainsQUAD + nDomainsQUADQUAD] + xMin[i + nDomainsQUAD + nDomainsQUADQUAD]+ngcdX[i + nDomainsQUAD + nDomainsQUADQUAD]))< eps0) &&
			(std::abs(xCoord[k][1]- (jjdY[i + nDomainsQUAD + nDomainsQUADQUAD] + yMin[i + nDomainsQUAD + nDomainsQUADQUAD]+ngcdY[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]- (iidX[i] + xMin[i]+ngcdX[i]))< eps0) &&
			(std::abs(xCoord[k][1]- (jjdY[i] + yMin[i]+ngcdY[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(2ngc + 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(2ngc + 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(2ngc + 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(2ngc + 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(2ngc + 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(2ngc + 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(2ngc + 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(2ngc + 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(2ngc + 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[2j] = xCoord[cells[jj][j]][0] + ngcdX[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(2ngc + 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(2ngc + 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[2j] = xCoord[cells[jj][j]][0] - ngcdX[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(2ngc + 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(2ngc + 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[2j + 1] = xCoord[cells[jj][j]][1] + ngcdY[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(2ngc + 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(2ngc + 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[2j + 1] = xCoord[cells[jj][j]][1] - ngcdY[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(2ngc + 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]-ngcdX[iQUAD])) < eps0&&std::fabs(yMin_temp - (yMax[iQUAD]-ngcdX[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]+ngcdX[iQUAD]))< eps0&&std::fabs(yMin_temp - (yMax[iQUAD]-ngcdX[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]+ngcdX[iQUAD]))< eps0&&std::fabs(yMax_temp - (yMin[iQUAD]+ngcdX[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(2ngc + 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] + 2n);
			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[2i + nDomainsTRIQUAD][3j + 0] + submeshCell0[i0 + QUADQUAD2QUAD[i][0]]);
			bcMap.push_back(bcMap_loc[2i + nDomainsTRIQUAD][3j + 1] + submeshCell0[i1 + i]);
			bcMap.push_back(-1);
			}

			// QUAD 2 QUADQUAD
			for (int j = 0; j < n; j++){
			bcMap.push_back(bcMap_loc[2i + nDomainsTRIQUAD][3(j+n) + 0] + submeshCell0[i1 + i]);
			bcMap.push_back(bcMap_loc[2i + 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[2i + 1 + nDomainsTRIQUAD][3j + 0] + submeshCell0[i0 + QUADQUAD2QUAD[i][1]]);
			bcMap.push_back(bcMap_loc[2i + 1 + nDomainsTRIQUAD][3j + 1] + submeshCell0[i1 + i]);
			bcMap.push_back(-1);
			}

			// QUAD 2 QUADQUAD
			for (int j = 0; j < n; j++){
			bcMap.push_back(bcMap_loc[2i + 1 + nDomainsTRIQUAD][3(j+n) + 0] + submeshCell0[i1 + i]);
			bcMap.push_back(bcMap_loc[2i + 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 + 2nDomainsQUADQUAD][3j + 0] + submeshCell0[i0 + RT2QUADQUAD[i][l]]);
			bcMap2.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2nDomainsQUADQUAD][3j + 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 + 2nDomainsQUADQUAD][3(j+n) + 0] + submeshCell0[i1 + i] + n);
			bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2nDomainsQUADQUAD][3(j+n) + 1] + submeshCell0[i0 + RT2QUADQUAD[i][l]]);
			bcMap.push_back(-1);
			bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2nDomainsQUADQUAD][3(j+n) + 0] + submeshCell0[i1 + i] + n + n0);
			bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2nDomainsQUADQUAD][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 + 2nDomainsQUADQUAD][3(j) + 0] + submeshCell0[i1 + i] + n);
			bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2nDomainsQUADQUAD][3(j) + 1] + submeshCell0[i0 + RT2QUAD[i][l]]);
			bcMap.push_back(-1);
			bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2nDomainsQUADQUAD][3(j) + 0] + submeshCell0[i1 + i] + n + n0);
			bcMap.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2nDomainsQUADQUAD][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 + 2nDomainsQUADQUAD][3j + 0] + submeshCell0[i0 + RT2QUADQUAD[i][l]]);
			bcMap2.push_back(bcMap_loc[indRT + nDomainsTRIQUAD + 2nDomainsQUADQUAD][3j + 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 (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

ConSol/ConSol/ConSol/Mesh/FastHybridMesh.cpp031392caa76amaster

FastHybridMesh.cppView Options

ConSol/ConSol/ConSol/Mesh/FastHybridMesh.cpp
031392caa76amaster

FastHybridMesh.cpp
View Options