diff --git a/ConSol/ConSol/Reconstr/HybridRBFWENO.cpp b/ConSol/ConSol/Reconstr/HybridRBFWENO.cpp
index f3dc07c..f8dd1d8 100644
--- a/ConSol/ConSol/Reconstr/HybridRBFWENO.cpp
+++ b/ConSol/ConSol/Reconstr/HybridRBFWENO.cpp
@@ -1,287 +1,287 @@
 //
 // Created by Fabian Moenkeberg on 2019-10-21.
 //
 
 #include "HybridRBFWENO.hpp"
 #include "ENO_2D_otf.hpp"
 #include "MRWENO_c.hpp"
 #include "CWENO_2D.h"
 #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"
 #include "Utility/Hesthaven/stdWenoWeights.hpp"
 
 void HybridRBFWENO::initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Source* source, Configuration config){
 
     if (config.getStencilChoice() == SYMMc){
         this->ENO = new class MRWENO_c();
     }else if (config.getStencilChoice() == SYMM){
         this->ENO = new class CWENO_2D();
     }else{
         this->ENO = new class ENO_2D_otf();
     }
     ENO->initialize(model, interpBase, limiter, mesh, source, config);
     this->model = model;
     this->source = source;
     this->limiter = limiter;
     this->interpBase = interpBase;
     this->order = ceil((config.getOrder()+1)/2.0)+1;
     this->nThreads = config.getNthreads();
     int n = config.getOrder();
     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 (config.getOrder() < 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";
 
     xCoord = mesh->getXCoord();
     ngc = mesh->getNgc();
 //    nX = xCoord[0].size()-1-2*ngc; // Number of cells in x-dir.
 //    nY = xCoord[1].size()-1-2*ngc; // Number of cells in y-dir.
 
     // Initialize linear weights WENO reconstruction to get higher order by combining the low
     // order reconstructions v_{i+1/2} = sum_{r=0}^{k-1} dw(r)v^{r}_{i+1/2}.
     dw = LinearWeights(order,0);
     Crec.resize(order + 1, std::vector<double> (order + 1));
     for (int r = -1; r < order; r++){
         Crec[r+1] = ReconstructWeights(order, r);
     }
 
     std::vector<std::vector<double>> beta0;
     beta.resize(order, std::vector<std::vector<double>>(order, std::vector<double>(order,0)));
     std::vector<double> xL;
     for (int r = 0; r < order; r++){
         xL.resize(order + 1);
         for (int ii = 0; ii <= order; ii++){
             xL[ii] = -0.5 + ii - r;
         }
 
         beta0 = betarcalc(xL, order);
 
         for (int ii = 0; ii < order; ii ++){
             for (int jj = 0; jj < order; jj ++){
                 beta[ii][jj][r] = beta0[ii][jj];
             }
         }
     }
 
     nCells = mesh->getNcells();
     nQuads = QWeights.size();
     nDomainsQUAD = mesh->getNdomainsQUAD() + mesh->getNdomainsQUADQUAD() + mesh->getNdomainsRQ();
 //    nDomainsQUADQUAD = mesh->getNdomainsQUADQUAD();
 //    nDomainsRQ = mesh->getNdomainsRQ();
     nDomainsTRIQUAD = mesh->getNdomainsTRIQUAD();
     nDomainsTRI = mesh->getNdomainsTRI();
     nDomainsRT = mesh->getNdomainsRT();
     submeshCell0 = mesh->getSubmeshCell0();
     submeshEdg0 = mesh->getSubmeshEdg0();
     nX = mesh->getNx();
     nY = mesh->getNy();
 
     // Quadrilateral weights
     order_scheme = config.getOrder();
     w0.resize(n/2 + (n%2));
     abs.resize(n/2 + (n%2));
     x_legendre.resize(n);
     weights_legendre.resize(n);
     w = new double[n];
     x = new double[n];
 
     legendre_set(n,x,w);
     if (order_scheme < 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] = w[i];
         }
     }
 
     delete [] w;
     delete [] x;
 }
 
 void HybridRBFWENO::reconstruct(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, std::vector<std::vector<double>> &U,
         MeshBase *mesh, double t, std::vector<std::vector<double>> &Source){
     int dimSyst = mesh->getNelem();
 
 
     ENO->reconstruct(UReconstr, U, mesh, t, Source);
 
     for (int ii = 0; ii < nDomainsQUAD; ii++){
     omp_set_num_threads(nThreads);
 #pragma omp parallel for
         for (int i = 0; i < nX[ii]; i++) {
             std::vector<std::vector<double>> uRec(dimSyst, std::vector<double>(2));
             std::vector<double> uloc(2*(order-1) + 1);
             std::vector<double> rec(2);
 
             for (int j = 0; j < nY[ii]-1; j++) {
                 for (int k = 0; k < dimSyst; k++) {
                     for (int l = 1 - order; l < order; l++){
                         uloc[l + order - 1] = U[i+ngc+(l+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k];
                     }
                     WENO_ev(uloc, rec);
                     uRec[k] = rec;
-                    limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
 //                    UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
 //                    UReconstr[1][i+(j+1)*nX[ii]  + submeshEdg0[ii + 1]][0][k] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
                 }
+                limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
                 for (int k = 0; k < dimSyst; k++) {
                     UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
                     UReconstr[1][i+(j+1)*nX[ii]  + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
                 }
             }
              int j = nY[ii]-1;
 
             for (int k = 0; k < dimSyst; k++) {
                 for (int l = 1 - order; l < order; l++){
                     uloc[l + order - 1] = U[i+ngc+(l+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k];
                 }
                 WENO_ev(uloc, rec);
                 uRec[k] = rec;
-                limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
 //                UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
 //                UReconstr[0][i+(j+1)*nX[ii]  + submeshEdg0[ii + 1]][0][k] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
             }
+            limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
             for (int k = 0; k < dimSyst; k++) {
                 UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
                 UReconstr[1][i+(j+1)*nX[ii]  + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
             }
         }
         omp_set_num_threads(nThreads);
 #pragma omp parallel for
         for (int j = 0; j < nY[ii]; j++) {
             std::vector<std::vector<double>> uRec(dimSyst, std::vector<double>(2));
             std::vector<double> uloc(2*(order-1) + 1);
             std::vector<double> rec(2);
 
             for (int i = 0; i < nX[ii]-1; i++) {
 
                 for (int k = 0; k < dimSyst; k++) {
                     for (int l = 1 - order; l < order; l++){
                         uloc[l+order-1] = U[i + ngc + l + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k];
                     }
                     WENO_ev(uloc, rec);
                     uRec[k] = rec;
-                    limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
 //                    UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
 //                    UReconstr[1][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
                 }
+                limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
                 for (int k = 0; k < dimSyst; k++) {
                     UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
                     UReconstr[1][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
                 }
             }
 
             int i = nX[ii] - 1;
             for (int k = 0; k < dimSyst; k++) {
                 for (int l = 1 - order; l < order; l++){
                     uloc[l+order-1] = U[i + ngc + l + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]][k];
                 }
                 WENO_ev(uloc, rec);
                 uRec[k] = rec;
-                limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
 //                UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
 //                UReconstr[0][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k];
             }
+            limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]);
             for (int k = 0; k < dimSyst; k++) {
                 UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][0];
                 UReconstr[1][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][1];
             }
         }
     }
 
     mesh->updateEdges(UReconstr, t);
 
 }
 
 std::vector<int> HybridRBFWENO::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 HybridRBFWENO::WENO_ev(std::vector<double> uloc, std::vector<double> &recloc){
 
     double vareps = 1e-6;
 //    std::vector<double> upl(order,0);
 //    std::vector<double> uml(order,0);
 //    std::vector<double> betar00(order,0);
     double upl;
     double uml;
     double betar;
     double tmp;
     recloc[0] = 0;
     recloc[1] = 0;
 //    std::vector<double> alphap(order);
 //    std::vector<double> alpham(order);
 
     double sumalpham = 0;
     double sumalphap = 0;
 
     if (order == 1){
         recloc[0] = uloc[0];
         recloc[1] = uloc[0];
     }else{
 
 
 //        std::vector<double> umh(order);
 
         for (int r = 0; r < order; r++){
             betar = 0;
             upl = 0;
             uml = 0;
             for (int j = 0; j < order; j++){
                 upl += Crec[r+1][j]*uloc[order-1 - r + j];
                 uml += Crec[r][j]*uloc[order-1 - r + j];
                 for (int i = 0; i < order; i++){
                     betar += uloc[order-1 - r + j]*beta[j][i][r]*uloc[order-1 - r + i];
                 }
 
             }
             tmp = (vareps+betar);
             sumalpham += dw[order-1-r]/(tmp*tmp);
             recloc[0]+= dw[order-1-r]/(tmp*tmp)*uml;
             sumalphap += dw[r]/(tmp*tmp);
             recloc[1]+= dw[r]/(tmp*tmp)*upl;
         }
 
         recloc[0] = recloc[0]/sumalpham;
         recloc[1] = recloc[1]/sumalphap;
 
 //        recloc[0] = uloc[1];
 //        recloc[1] = uloc[1];
     }
 }
 
 void HybridRBFWENO::exportStencils(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, std::string output_file){
     ENO->exportStencils(U, mesh, t, output_file);
 }
\ No newline at end of file