diff --git a/ConSol/ConSol/Reconstr/ENO_2D.cpp b/ConSol/ConSol/Reconstr/ENO_2D.cpp
index a4c9145..6181a63 100644
--- a/ConSol/ConSol/Reconstr/ENO_2D.cpp
+++ b/ConSol/ConSol/Reconstr/ENO_2D.cpp
@@ -1,375 +1,375 @@
 //
 // Created by Fabian Moenkeberg on 30.07.18.
 //
 
 #include "ENO_2D.hpp"
 #include <assert.h>
 #include "../Utility/MatrixVectorProd.hpp"
 //#include "omp.h"
 #include <omp.h>
 #include <numeric>
 #include "../Utility/ddGeneralizedNewton.hpp"
 #include "../Utility/signPropFulfilled.hpp"
 #include "../Utility/isWellDef.hpp"
 #include "Utility/gauss-legendre-quadrature/gauss_legendre.hpp"
 #include "H5Cpp.h"
 
 void ENO_2D::initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Configuration config){
     this->model = model;
     this->limiter = limiter;
     this->interpBase = interpBase;
     this->order = config.getOrder();
     this->nThreads = config.getNthreads();
 	int n = order;
     std::vector<double> w0(n/2 + (n%2));
     std::vector<double> abs(n/2 + (n%2));
     Qlambdas.resize(n, std::vector<double>(2));
     QWeights.resize(n);
     double *w = new double[n];
     double *x = new double[n];
 
     legendre_set(n,x,w);
     if (order < 2){
         Qlambdas[0][0] = 0.5; Qlambdas[0][1]= 0.5;
         QWeights[0] = 1;
     }else {
         for (int i = 0; i < n; i++) {
             Qlambdas[i][0] = (1 + x[i]) / 2;
             Qlambdas[i][1] = (1 - x[i]) / 2;
             QWeights[i] = w[i] / 2;
         }
     }
 
     if (nThreads == 0){
         nThreads = omp_get_max_threads();
     }
 
     std::cout << "  Reconstruction: Running threads: " << nThreads << " out of " << omp_get_max_threads() << "\n";
 
     switch (config.getDegPoly()+1)
     {
         case 2:
             nStencil = 5;
             maxLevelNeighbour = 3;
 //            nStencil = 7;
 //            maxLevelNeighbour = 2;
             break;
         case 3:
-            nStencil = 14;
-            maxLevelNeighbour = 6;
+            nStencil = 12;
+            maxLevelNeighbour = 5;
 //            nStencil = 14;
 //            maxLevelNeighbour = 3;
             break;
         case 4:
             nStencil = 23;
             maxLevelNeighbour = 8;
             break;
         default:
             nStencil = config.getOrder()*(config.getOrder()+1)/2 + 1;
             maxLevelNeighbour = ceil(double(nStencil)/2);
             break;
     }
 
     delete [] w;
     delete [] x;
 }
 
 void ENO_2D::exportStencils(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, std::string output_file){
     int dimSyst = mesh->getNelem();
     int nCells = mesh->getNcells();
 
     std::vector<std::vector<int>> saveStencil(nCells,std::vector<int>(nStencil));
 
     omp_set_num_threads(nThreads);
 #pragma omp parallel for
     for (int i = 0; i < nCells; i++) {
 
         int indInternal = mesh->getInternal(i);
 
         for (int k = 0; k < dimSyst; k++) {
             std::vector<int> stencil = this->getStencil(indInternal, mesh, &U, k);
             saveStencil[i] = stencil;
         }
     }
     exportVec(saveStencil, mesh, U, output_file);
 }
 
 
 void ENO_2D::reconstruct(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, std::vector<std::vector<double>> &U, MeshBase *mesh, double t){
     int dimSyst = mesh->getNelem();
     int nCells = mesh->getNcells();
     int nQuads = QWeights.size();
 
     std::vector<std::vector<int>> saveStencil(nCells,std::vector<int>(nStencil));
     omp_set_num_threads(nThreads);
 #pragma omp parallel for
     for (int i = 0; i < nCells; i++) {
 
         int indInternal = mesh->getInternal(i);
 //        if (indInternal==105){
 //            std::cout << i << "\n";
 //        }
         std::vector<int> c2edg = mesh->getC2E(indInternal);
         std::vector<std::vector<double>> evPoints(NODE_NUM * nQuads, std::vector<double>(2, 0));
         std::vector<std::vector<double>> edg_i(2,std::vector<double>(2));
 
         for (int l = 0; l < NODE_NUM; l++) {
             for (int j = 0; j < nQuads; j++) {
                 edg_i = mesh->edge2Vertex(c2edg[l]);
                 evPoints[NODE_NUM * j + l][0] = Qlambdas[j][0] * edg_i[0][0] + Qlambdas[j][1] * edg_i[0][1]; // x-Coordinate
                 evPoints[NODE_NUM * j + l][1] = Qlambdas[j][0] * edg_i[1][0] + Qlambdas[j][1] * edg_i[1][1]; // y-Coordinate
             }
          }
 
         std::vector<std::vector<double>> uRec(dimSyst, std::vector<double>(NODE_NUM * nQuads));
         for (int k = 0; k < dimSyst; k++) {
             std::vector<int> stencil = this->getStencil(indInternal, mesh, &U, k);
 //            std::vector<int> stencil{1,2,3,4,5};
             saveStencil[i] = stencil;
             int sizeStencil = (int) stencil.size();
             std::vector<double> locValues(sizeStencil);
             for (int ii = 0; ii < sizeStencil; ii++) {
                 locValues[ii] = U[stencil[ii]][k];
             }
 			std::vector<double> intValues;
 			if (stencil.size() <= 1){
 				int nEv = evPoints.size();
 				intValues.resize(nEv);
 				for (int l = 0; l < nEv; l++) {
                     	intValues[l] = locValues[0];
          		}
 			}else{
             	intValues = interpBase->evaluate(evPoints, mesh, stencil, locValues);
 			}
             uRec[k] = intValues;
 //            for (int jj = 0; jj < intValues.size(); jj++){
 //                uRec[k][jj] = U[indInternal][k];
 //            }
         }
         for (int j = 0; j < NODE_NUM; j++) {
             if (mesh->getRight(c2edg[j]) == indInternal) {
                 for (int l = 0; l < nQuads; l++) {
                     for (int k = 0; k < dimSyst; k++) {
                         UReconstr[0][c2edg[j]][l][k] = uRec[k][NODE_NUM * l + j];
                     }
                 }
 
             } else {
                 for (int l = 0; l < nQuads; l++) {
                     for (int k = 0; k < dimSyst; k++) {
                         UReconstr[1][c2edg[j]][l][k] = uRec[k][NODE_NUM * l + j];
                     }
                 }
             }
         }
     }
 
     mesh->updateEdges(UReconstr, t);
 }
 
 std::vector<int> ENO_2D::getStencil(int cellIndex, MeshBase *mesh, std::vector<std::vector<double>> *U, int k){
     int indMinSmInd = 0;
     std::vector<std::vector<int>> stencils = mesh->getStencils(cellIndex);
     return stencils[indMinSmInd];
 }
 
 
 void ENO_2D::exportVec(std::vector<std::vector<int>> vec, MeshBase *mesh, std::vector<std::vector<double>> U, std::string output_path){
     std::string output_file = output_path + "/stencils.h5";
     H5::H5File *file = new H5::H5File(output_file, H5F_ACC_TRUNC);
 
     hsize_t dimsU[2];
     int nSpace = (int) U.size();
     int nElem = (int) U[0].size();
 
 
     dimsU[0] = nSpace;
     dimsU[1] = 1;
 
     H5::DataSpace dataspaceU(2,dimsU);
 
 
     H5::DataType datatypeU(H5::PredType::NATIVE_DOUBLE);
     double* U1Vec0 = new double [nSpace];
     for (int i = 0; i < nElem; i++){
 
         char name[10];
         sprintf(name, "variable %d",i);
         H5::DataSet dataset = file->createDataSet(name, datatypeU, dataspaceU);
 
         for (int k = 0; k < nSpace; k++){
             U1Vec0[k] = U[k][i];
         }
         for (int i = 0; i < 1; i++){
             dataset.write(U1Vec0,H5::PredType::NATIVE_DOUBLE);
         }
 
         dataset.close();
     }
 
 
     hsize_t dimsf[2];
     int nVec1= (int) vec.size();
     int nVec2 = (int) vec[0].size();
 
 
     dimsf[0] = nVec1;
     dimsf[1] = nVec2;
 
     H5::DataSpace dataspace(2,dimsf);
 
     H5::DataType datatype(H5::PredType::NATIVE_INT);
     int* U1Vec = new int [nVec1*nVec2];
     H5::DataSet dataset = file->createDataSet("vector", datatype, dataspace);
 
     for (int j = 0; j < nVec1; j++){
         for (int k = 0; k < nVec2; k++){
             U1Vec[k+j*nVec2] = vec[j][k];
         }
     }
 
     dataset.write(U1Vec,H5::PredType::NATIVE_INT);
     dataset.close();
     dataspace.close();
 
     hsize_t dimsXC[2];
     H5::Group Mesh = file->createGroup("mesh");
     std::vector<std::vector<double>> xCoord = mesh->getXCoord();
 //    int nXCoord = (int) xCoord.size();
     int dim = mesh->getNdims();
     std::vector<std::vector<int>> cells = mesh->getCells();
 
     if (mesh->getMeshType() == HYBRID){
         std::vector<double> dX = mesh->getDx();
         std::vector<double> dY = mesh->getDy();
         std::vector<double> xMin = mesh->getXmin();
         std::vector<double> yMin = mesh->getYmin();
         std::vector<int> Nx = mesh->getNx();
         std::vector<int> Ny = mesh->getNy();
         int ngc = mesh->getNgc();
         int nXCoord2;
         int nYCoord2;
         std::vector<std::vector<double>> xCoord2;
         std::vector<std::vector<int>> cells2;
         for (int ii = 0; ii < dX.size(); ii++){
             nXCoord2 = Nx[ii]+1;
             nYCoord2 = Ny[ii]+1;
 //            std::cout << "test" << ", ";
             xCoord2.resize((nXCoord2 + 2*ngc) * (nYCoord2 + 2*ngc), std::vector<double>(2));
 
             for (int i = 0; i < nXCoord2 + 2*ngc; i++) {
                 for (int j = 0; j < nYCoord2 + 2*ngc; j++) {
                     xCoord2[i + j * (nXCoord2 + 2*ngc)][0] = xMin[ii] + i*dX[ii];
                 }
             }
 
             for (int i = 0; i < nXCoord2 + 2*ngc; i++) {
                 for (int j = 0; j < nYCoord2 + 2*ngc; j++) {
                     xCoord2[i + j * (nXCoord2 + 2*ngc)][1] = yMin[ii] + j*dY[ii];
                 }
             }
 
             int n1 = (int) (nXCoord2-1 + 2*ngc)*(nYCoord2-1+2*ngc);
             int n2 = (int) 4;
 //        std::vector<std::vector<int>> cells();
             cells2.resize(n1, std::vector<int>(n2));
 
             for (int ix = 0; ix < nXCoord2 + 2*ngc - 1; ix++) {
                 for (int iy = 0; iy < nYCoord2 + 2*ngc - 1; iy++) {
                     cells2[(ix + iy*(nXCoord2 - 1 + 2*ngc))][0] = ix + iy*(nXCoord2 + 2*ngc) + xCoord.size();
                     cells2[(ix + iy*(nXCoord2 - 1 + 2*ngc))][1] = ix+1 + iy*(nXCoord2 + 2*ngc) + xCoord.size();
                     cells2[(ix + iy*(nXCoord2 - 1 + 2*ngc))][2] = ix+1 + (iy+1)*(nXCoord2 + 2*ngc) + xCoord.size();
                     cells2[(ix + iy*(nXCoord2 - 1 + 2*ngc))][3] = ix + (iy+1)*(nXCoord2 + 2*ngc) + xCoord.size();
                 }
             }
 
             cells.insert(cells.end(), cells2.begin(), cells2.end());
             xCoord.insert(xCoord.end(), xCoord2.begin(), xCoord2.end());
         }
     }
     int nXCoord = (int) xCoord.size();
 
     dimsXC[0] = dim;
     dimsXC[1] = nXCoord;
 
     double* xCoordVector = new double[nXCoord * dim];
 
     for (int i = 0; i < nXCoord; i++){
         for (int j = 0; j < dim; j++){
             xCoordVector[i + j * nXCoord] = xCoord[i][j];
         }
     }
 
     H5::DataSpace dataspaceXCoord(2,dimsXC);
     H5::DataType datatypeXCoord(H5::PredType::NATIVE_DOUBLE);
 
     H5::DataSet dataSetXCoord = Mesh.createDataSet("VCoordinates", datatypeXCoord, dataspaceXCoord);
     dataSetXCoord.write(xCoordVector, H5::PredType::NATIVE_DOUBLE);
     dataSetXCoord.close();
     dataspaceXCoord.close();
 
     int n1 = (int) cells.size();
     int n2 = (int) cells[0].size();
     hsize_t dimsC[2];
     dimsC[0] = n1;
     dimsC[1] = 4;
 
     H5::DataSpace dataspaceCells(2,dimsC);
     H5::IntType datatypeCells(H5::PredType::NATIVE_INT);
 
     int * cellsVector = new int[n1*4];
     std::vector<int> cellsVector0(0);
 //    std::cout << n1 << "," << n2 << "\n";
     for (int i = 0; i < n1; i++){
         std::vector<int> cell_tmp = cells[i];
 //        for (int j = 0; j < n2;j++){
 //            std::cout << cells[i][j]<< ", ";
 //            cellsVector[i+j*n1] = cell_tmp[j];
             if (cell_tmp.size()==3){
                 cell_tmp.push_back(cell_tmp[2]);
             }
             cellsVector0.insert(cellsVector0.end(),cell_tmp.begin(),cell_tmp.end());
 //            std::cout << cell_tmp[j]<< ", ";
 //        }
 //        std::cout << "\n";
     }
 //    std::cout << cellsVector[0]<< ","<< cellsVector[1]<< ","<< cellsVector[2]<< ","  << "\n";
 
     H5::DataSet dataSetCells = Mesh.createDataSet("Cells", datatypeCells, dataspaceCells);
     dataSetCells.write(cellsVector0.data(), H5::PredType::NATIVE_INT);
     dataSetCells.close();
     dataspaceCells.close();
 
     std::vector<int> internal = mesh->getInternal();
     int* internalVector = internal.data();
 
     hsize_t dimsI[2];
     dimsI[0] = 1;
     dimsI[1] = internal.size();
 
     H5::DataSpace dataspaceI(2,dimsI);
     H5::DataType datatypeI(H5::PredType::NATIVE_INT);
 
     H5::DataSet dataSetI = Mesh.createDataSet("InternalCells", datatypeI, dataspaceI);
     dataSetI.write(internalVector, H5::PredType::NATIVE_INT);
     dataSetI.close();
     dataspaceI.close();
     H5::DataSpace dspace(H5S_SCALAR);
     H5::Attribute attMesh = Mesh.createAttribute("nDimension", H5::PredType::NATIVE_INT,dspace);
     int ndims = mesh->getNdims();
     attMesh.write(H5::PredType::NATIVE_INT,&ndims);
     attMesh.close();
 
     attMesh = Mesh.createAttribute("dimSystem", H5::PredType::NATIVE_INT,dspace);
     int nelem = mesh->getNelem();
     attMesh.write(H5::PredType::NATIVE_INT,&nelem);
     attMesh.close();
     Mesh.close();
 
 //
     file->close();
     delete [] U1Vec;
     delete [] xCoordVector;
     delete [] cellsVector;
     delete file;
 }