diff --git a/ConSol/ConSol/Mesh/Tri.cpp b/ConSol/ConSol/Mesh/Tri.cpp
index 5487e26..dfe2a8b 100644
--- a/ConSol/ConSol/Mesh/Tri.cpp
+++ b/ConSol/ConSol/Mesh/Tri.cpp
@@ -1,639 +1,639 @@
 //
 // Created by Fabian Moenkeberg on 18.07.18.
 //
 
 #include "Tri.h"
 #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>
 
 # define NODE_NUM 3
 
 using namespace std;
 
 void Tri::initialize(BcBase *bc, Configuration config) {
     this->bc = bc;
     int internal_num;
     int bcMap_num;
     int edg_num;
     int element_order;
     nDomainsQUAD = 0;
     nDomainsTRI = 1;
     nDomainsTRIQUAD = 0;
     bcMap2.resize(0);
 
     string gmsh_filename = config.getGridFile();
     int node_num;
     this->nThreads = config.getNthreads();
 
     if (nThreads == 0){
         nThreads = omp_get_max_threads();
     }
 
     std::cout << "  Triangulation: Running threads: " << nThreads << " out of " << omp_get_max_threads() << "\n";
 
     cout << "\n";
     cout << "  Load mesh data from a file.\n";
 
-    if (false){
+    if (true){
     //
     //  Get the data size.
     //
     //    gmsh_io::gmsh_size_read( gmsh_filename, node_num, ndims, nCells,
     //                             element_order );
         gmsh_io::mesh_size_read( gmsh_filename, node_num, ndims, nCells,
                                  element_order, internal_num , bcMap_num, edg_num);
     //
     //  Print the sizes.
     //
         cout << "\n";
         cout << "  Node data read from file \"" << gmsh_filename << "\"\n";
         cout << "\n";
         cout << "  Number of nodes = " << node_num << "\n";
         cout << "  Spatial dimension = " << ndims << "\n";
         cout << "  Number of elements = " << nCells << "\n";
         cout << "  Element order = " << element_order << "\n";
 
     //
     //  Get the data.
     //
         gmsh_io::mesh_data_read ( gmsh_filename, ndims, node_num, xCoord,
                                   3, nCells, cells , internal, bcMap, e2c, edg, stencils,
                                   neighbours, c2e, edgesAtBC);
 
         nCells = internal.size();
         nCellsTot = cells.size();
         nedges = edg.size();
 
 
         for (int i = 0; i < edgesAtBC.size(); i++){
             edgesAtBC[i][2]++;
             edgesAtBC[i].push_back(edgesAtBC[i][1]);
             if (config.getBc() == PERIODIC){
                 edgesAtBC[i][1] = 4;
             }else{
                 assert(false && "TRI: not implemented BC!");
             }
 
         }
         std::vector<int> bcMap_tmp(0);
         for (int i = 0; i < bcMap.size(); i++){
             if(bcMap[i]!=i){
                 bcMap_tmp.push_back(i);
                 bcMap_tmp.push_back(bcMap[i]);
                 bcMap_tmp.push_back(-1);
             }
         }
 
-
+        bcMap = bcMap_tmp;
         bcMap.resize(0);
     }else{
         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();
         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::mesh::getElementsByType(eleType2D, elementTags, nodeTags, s);
                 gmsh::logger::write("- " + std::to_string(elementTags.size()) +
                                     " elements in surface " + std::to_string(s));
 
                 // 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);
 
         int nNodes = coord.size()/3;
         xCoord.resize(nNodes,std::vector<double>(2));
 
         for (int i = 0; i < nNodes; i++){
             xCoord[i][0] = coord[3*i];
             xCoord[i][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 == "absorb"){
                 listPhysicalNameEdges.push_back(ABSORBINGBC);
             }else{
                 assert(false && "Read Mesh: 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<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 == "tri") {
                 listPhysicalNameElements.push_back(TRI);
             }else if (physicalName == "quad"){
                 listPhysicalNameElements.push_back(QUAD);
             }else{
                 assert(false && "Read Mesh: 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. ");
         }
 
         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;
         for (int i = 0; i < nElementTypes; i++) {
             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);
             int nElements = elementTags[i].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);
                 c2e.push_back(cell0);
                 cells.push_back(cell0);
                 for (int j = 0; j <numNodes; j++){
                     c2e[iElement][j]  = -1;
                     cells[iElement][j] = nodeTagsElements[i][numNodes*ii  + j]-1;
                     face[0] =nodeTagsElements[i][numNodes*ii  + map[2*j]]-1;
                     face[1] =nodeTagsElements[i][numNodes*ii  + map[2*j + 1]]-1;
 
                     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;
                         e2c[indFace][1] = iElement;
                     } else{
                         nFaces++;
                         c2e[iElement][j] = nFaces-1;
                         edg.push_back(face);
                         e2c.push_back({iElement,-1});
                     }
                     faceExists = 0;
 
                 }
             }
             int nF = edg.size();
             int nF2c = e2c.size();
             int nc2e = c2e.size();
 
             nElements0 += numNodes*nElements;
         }
 
         // Calculate neighbours mapping or neighbour mapping
         nElements0 = 0;
         iElement = -1;
         for (int i = 0; i < nElementTypes; i++) {
             int n;
             int nElements = elementNumbers[i];
             int numNodes = elementSizes[i];
             std::vector<int> cell0(numNodes);
 
             for (int ii = 0; ii < nElements; ii++){
                 iElement++;
                 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;
                     }
                 }
             }
 
             nElements0 += numNodes*nElements;
         }
 
         // 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 = (coord[3*cells[i][(j+1)%numNodes]] - coord[3*cells[i][j]])*(coord[3*cells[i][(j+2)%numNodes]+1] - coord[3*cells[i][(j+1)%numNodes]+1]) -
                      (coord[3*cells[i][(j+1)%numNodes]+1] - coord[3*cells[i][j]+1])*(coord[3*cells[i][(j+2)%numNodes]] - coord[3*cells[i][(j+1)%numNodes]]);
                 if (cp < 0){
                     assert(false && "Grid Generation: 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;
             }
         }
 
         // Create "right" cells of each edge and create edgeAtBCMap
         // Note: for each edge "left" cell is always in domain.
         edgesAtBC.resize(0, std::vector<int>(0));
         std::vector<int> edgeInfo(3);
         int nBC = listPhysicalNodeTagsEdges.size();
         for (int i = 0; i < nFaces; i++){
             if (e2c[i][1] == -1){
                 edgeInfo[0] = i;
                 for (int j = 0; j < nBC; j++){
                     int n = listPhysicalNodeTagsEdges[j].size();
                     int count = 0;
                     for (int ii = 0; ii < n; 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;
                 edgesAtBC.push_back(edgeInfo);
             }
         }
 
 
 
         nCellsTot = cells.size();
         nCells = nCellsTot;
         nedges = edg.size();
 
         internal.resize(nCells);
         bcMap.resize(0);
         for (int i = 0; i < nCells; i++){
             internal[i] = i;
             internal[i] = i;
         }
     }
 
     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 && "Tri: model not defined");
             break;
     };
 
     normals.resize(nedges,std::vector<double>(ndims));
     std::vector<double> tmp(ndims);
     double m;
     for (int i = 0; i < nedges; i++){
         tmp[0] = xCoord[edg[i][1]][0]-xCoord[edg[i][0]][0];
         tmp[1] = xCoord[edg[i][1]][1]-xCoord[edg[i][0]][1];
         m = std::sqrt(tmp[0]*tmp[0]+tmp[1]*tmp[1]);
         normals[i][0] = tmp[1]/m;
         normals[i][1] = -tmp[0]/m;
     }
 
 
     size.resize(nCellsTot);
     dS.resize(nCellsTot,std::vector<double>(3));
 
     double s;
     std::vector<double> inRadius(nCellsTot);
 //    std::cout << "test" << nCellsTot << " \n";
     double alpha = 0;
     for (int i = 0; i < nCellsTot; i++){
 
         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);
         for (int j = 0; j < 3; j++) {
             tmp[0] = xCoord[cells[i][(j+1)%3]][0] - xCoord[cells[i][j]][0];
             tmp[1] = xCoord[cells[i][(j+1)%3]][1] - xCoord[cells[i][j]][1];
 //            std::cout << c2e[i][j] <<  "\n";
             dS[i][j] = (tmp[1] * normals[c2e[i][j]][0] - tmp[0] * normals[c2e[i][j]][1]) / size[i];
         }
 //        std::cout << i <<  "\n";
     }
 //    std::cout << "test \n";
     dxMin.resize(2);
     double minInradius = *std::min_element(std::begin(inRadius), std::end(inRadius));
     dxMin[0] = minInradius/alpha;
     dxMin[1] = minInradius/alpha;
 
 //    int nEdgesBC = indBCLeftIsOut.size()+indBCRightIsOut.size();
 //    edgesAtBC.resize(nEdgesBC,std::vector<int>(3));
 //    for (int i = 0; i < indBCLeftIsOut.size(); i++){
 //        edgesAtBC[i][0] = indBCLeftIsOut[i][0];
 //        edgesAtBC[i][1] = indBCLeftIsOut[i][1];
 //        edgesAtBC[i][2] = 1;
 //    }
 //    for (int i = 0; i < indBCRightIsOut.size(); i++){
 //        edgesAtBC[indBCLeftIsOut.size()+i][0] = indBCRightIsOut[i][0];
 //        edgesAtBC[indBCLeftIsOut.size()+i][1] = indBCRightIsOut[i][1];
 //        edgesAtBC[indBCLeftIsOut.size()+i][2] = 0;
 //    }
     submeshEdg0.resize(1,0);
     submeshEdg0.push_back(nedges);
     bc->initialize(this, config);
 }
 
 
 std::vector<double> Tri::MidPoint(int i){
     std::vector<double> midPts(ndims);
     midPts[0] = (xCoord[cells[i][0]][0] + xCoord[cells[i][1]][0] + xCoord[cells[i][2]][0])/3;
     midPts[1] = (xCoord[cells[i][0]][1] + xCoord[cells[i][1]][1] + xCoord[cells[i][2]][1])/3;
     return midPts;
 };
 
 void Tri::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){
 //    std::vector<std::vector<double>> ret(nCellsTot, std::vector<double>(nelem));
 //    std::vector<std::vector<double>> flux = numFluxFct->Feval(normals, UReconstr[1],UReconstr[0],t,dt);
 //    std::cout << flux.size() << "," << flux[0].size() << "\n";
 //    for (int i = 0; i < 20; i++){
 //        std::cout << flux[i][0]<< "\n";
 //    }
 
 //    std::vector<double> flux(nelem);
 //    int edge;
     omp_set_num_threads(nThreads);
 #pragma omp parallel for
     for (int i = 0; i < nCellsTot; i++){
 //        if (i == 6376){
 //            std::cout << i << "\n";
 //        }
         for (int k = 0; k < nelem; k++){
             numFlux[i][k] = 0;
         }
         for (int l = 0; l < 3; l++){
             int edge = c2e[i][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[i][k] -= flux[k]*dS[i][l];
             }
         }
     }
 };
 
 std::vector<std::vector<double> > Tri::evalFunc(AbstractFunct *Funct, evaluationMethod evalMethod){
 //    int nElem = Funct->getDim();
     std::vector<std::vector<double> > ret(nCellsTot, std::vector<double>(nelem,0));
 //    std::cout << "test \n";
     switch (evalMethod){
         case MEANVALUE:{
 
             double node_xy2[2*NODE_NUM];
             int i;
             int order_num;
             int rule;
             double *w;
             double *xy;
             double *xy2;
             rule = 8;
 
             order_num = dunavant_order_num ( rule );
             xy = new double[2*order_num];
             xy2 = new double[2*order_num];
             w = new double[order_num];
 
             dunavant_rule ( rule, order_num, xy, w );
             std::vector<double> x_ev(2);
             std::vector<double> F_ev;
 
             for (i = 0; i < nCellsTot; i++){
 
                 node_xy2[0] = xCoord[cells[i][0]][0];
                 node_xy2[1] = xCoord[cells[i][0]][1];
                 node_xy2[2] = xCoord[cells[i][1]][0];
                 node_xy2[3] = xCoord[cells[i][1]][1];
                 node_xy2[4] = xCoord[cells[i][2]][0];
                 node_xy2[5] = xCoord[cells[i][2]][1];
 
 
                 reference_to_physical_t3 ( node_xy2, order_num, xy, xy2 );
 
                 std::fill(ret[i].begin(), ret[i].end(), 0);
 //                std::cout << order_num << ", " << nElem << ", ";
                 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);
 //                    std::cout << F_ev.size() << "\n";
                     for (int l = 0; l<nelem; l++){
                         ret[i][l] += w[k]*F_ev[l];
                     }
                 }
 //                std::vector<double> m = MidPoint(i);
 //                ret[i] = Funct->F(MidPoint(i));
 //                std::cout << ret[i].size()<< ", " << nelem;
             }
             delete [] w;
             delete [] xy;
             delete [] xy2;
             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>> Tri::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;
 }
 # undef NODE_NUM
\ No newline at end of file