diff --git a/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.cpp b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.cpp
index 351870c..2b3f0c0 100644
--- a/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.cpp
+++ b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.cpp
@@ -1,329 +1,333 @@
 //
 // Created by moenkeberg on 11.02.20.
 //
 
 #include "Multiquadratics_2D_WENO_c.hpp"
 #include "cblas.h"
 
 Multiquadratics_2D_WENO_c::Multiquadratics_2D_WENO_c(int order, double eps){
     this->order = order;
     this->order_scheme = 2*((order+1)/2 + (order+1)%2);
     this->shapeParam = eps;
     this->dim = 2;
     this->degPoly = order -1;
     this->name = MULTIQUADS;
     setMulti_Index(degPoly+1);
 }
 
 Multiquadratics_2D_WENO_c::Multiquadratics_2D_WENO_c(int order, int degPoly, double eps){
     this->order = 1;
 //    this->order_scheme = 2*((order+1)/2 + (order+1)%2);
     this->order_scheme = order+1;
     this->shapeParam = eps;
     this->dim = 2;
     this->degPoly = degPoly;
     this->name = MULTIQUADS;
     setMulti_Index(degPoly+1);
 }
 
 std::vector<double> Multiquadratics_2D_WENO_c::evaluate(std::vector<std::vector<double>> evPoints, MeshBase *mesh, std::vector<std::vector<int>> stencil, std::vector<std::vector<double>> values){
     int nStencils = (int) stencil.size();
     int nStencil;
 
     double rad = 0;
     for (int i = 0; i < nStencils; i++){
         rad = std::max(stableRad / mesh->getApproxStencilSize(stencil[i]), rad);
     }
     std::vector<int> n = {1,4,7};
-    int degPolyMax = 2;
+
 
     double eps0 = 0.1/rad;
     double eps = rad;
     int nEv = (int) evPoints.size();
     std::vector<std::vector<double>> coeff(nStencils);
 
     std::vector<double> d;
     std::vector<std::vector<double>> xi(NODE_NUM,std::vector<double>(2));
     xi = mesh->getXCoord(mesh->getCells(stencil[0][0]));
     double C0 = calcSize(xi);
     double dx0 = (calcMidPt(xi,0));
     double dy0 = (calcMidPt(xi,1));
     double dx1;
     double dy1;
     double dx2;
     double dy2;
     int ns = 4;
     std::vector<double>a(3);
     std::vector<double>b(3);
     std::vector<double>beta(3);
     double a0;
     double b0;
     double smInd0 = 0;
     double tmp = 1;
     int itmp = 0;
+
     for (int ii = 1; ii < nStencils-1; ii++){
         nStencil = stencil[ii].size();
-        if (nStencil-1 == ns){
-            xi = mesh->getXCoord(mesh->getCells(stencil[ii][3]));
+        if (nStencil == ns){
+            xi = mesh->getXCoord(mesh->getCells(stencil[ii][2]));
             dx1 = (calcMidPt(xi,0));
             dy1 = (calcMidPt(xi,1));
-            xi = mesh->getXCoord(mesh->getCells(stencil[ii][4]));
+            xi = mesh->getXCoord(mesh->getCells(stencil[ii][3]));
             dx2 = (calcMidPt(xi,0));
             dy2 = (calcMidPt(xi,1));
-            a0 = (((values[ii][3] - values[ii][2])*(dy2 - dy0) - (values[ii][4] - values[ii][2])*(dy1 - dy0))/((dx2-dx0)*(dy1-dy0)-(dx1-dx0)*(dy2-dy0)));
-            b0 = (((values[ii][4] - values[ii][2])*(dy1 - dy0) - (values[ii][3] - values[ii][2])*(dy2 - dy0))/((dx2-dx0)*(dy1-dy0)-(dx1-dx0)*(dy2-dy0)));
+            a0 = (((values[ii][2] - values[ii][1])*(dy2 - dy0) - (values[ii][3] - values[ii][1])*(dy1 - dy0))/((dx2-dx0)*(dy1-dy0)-(dx1-dx0)*(dy2-dy0)));
+            b0 = (((values[ii][3] - values[ii][1])*(dy1 - dy0) - (values[ii][2] - values[ii][1])*(dy2 - dy0))/((dx2-dx0)*(dy1-dy0)-(dx1-dx0)*(dy2-dy0)));
             a[itmp] = a0;
             b[itmp] = b0;
             beta[itmp] = a0*a0 +b0*b0;
             itmp++;
         }
     }
+    int degPolyMax = 0;
 
     if (itmp == 3){
         tmp = (std::fabs(beta[0]-beta[1]) +  std::fabs(beta[0]-beta[2]) +  std::fabs(beta[2]-beta[1]))/3;
         double sigma0 = ((1 + tmp*tmp/(beta[0] + 1e-10))/3.0);
         double sigma1 = ((1 + tmp*tmp/(beta[1] + 1e-10))/3.0);
         double sigma2 = ((1 + tmp*tmp/(beta[2] + 1e-10))/3.0);
         smInd0 = my_pow(sigma0/(sigma0 +sigma1 +sigma2)*a[0],2) + my_pow(sigma0/(sigma0 +sigma1 +sigma2)*b[0],2)
                  + my_pow(sigma1/(sigma0 +sigma1 +sigma2)*a[1],2) + my_pow(sigma1/(sigma0 +sigma1 +sigma2)*b[1],2)
                  + my_pow(sigma2/(sigma0 +sigma1 +sigma2)*a[2],2) + my_pow(sigma2/(sigma0 +sigma1 +sigma2)*b[2],2);
+        degPolyMax = 2;
 
     }else if (itmp == 2){
         tmp = (std::fabs(beta[0]-beta[1]));
         double sigma0 = ((1 + tmp*tmp/(beta[0] + 1e-10))/2.0);
         double sigma1 = ((1 + tmp*tmp/(beta[1] + 1e-10))/2.0);
         smInd0 = my_pow(sigma0/(sigma0 +sigma1)*a[0],2) + my_pow(sigma0/(sigma0 +sigma1)*b[0],2)
                  + my_pow(sigma1/(sigma0 +sigma1)*a[1],2) + my_pow(sigma1/(sigma0 +sigma1)*b[1],2);
     }else if (itmp == 1){
         smInd0 = my_pow(a[0],2) + my_pow(b[0],2);
     }else{
         smInd0 = 1;
     }
 
     smInd0 *= C0;
+//    std::cout << smInd0 << std::endl;
 
     std::vector<double> smInd(3);
     std::vector<std::vector<double>> ret(3, std::vector<double>(nEv,0));
     std::vector<double> res(nEv,0);
     std::vector<double> coeff_old;
     int nRbf_old;
     int nP_old;
     int nRbf;
     std::vector<std::vector<std::vector<double>>> xB_;
     std::vector<double> values_;
     std::vector<double> values_old;
 
     std::vector<double> values_tmp;
 
     for (int ii = 0; ii <= degPolyMax; ii++) {
         nStencil = stencil[ii].size();
         std::vector<std::vector<std::vector<double>>> xB(nStencil, std::vector<std::vector<double>>(NODE_NUM,
                                                                                                     std::vector<double>(
                                                                                                             2)));
 
 
 
 //        values[2] = 14;
         if (ii == 0) {
             nRbf = 1;
             xB_.resize(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
             xi = mesh->getXCoord(mesh->getCells(stencil[0][0]));
             xB_[0] = xi;
             values_.resize(nRbf);
             values_[0] = values[0][0];
         }else if (ii == 1) {
             nRbf = 4;
             xB_.resize(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
             values_.resize(nRbf);
 
             for(int i = 0; i < nRbf; i++){
                 xi = mesh->getXCoord(mesh->getCells(stencil[nStencils-1][i]));
                 xB_[i] = xi;
                 values_[i] = values[nStencils-1][i];
             }
         }else{
             nRbf = 7;
 //            nRbf = 6;
             xB_.resize(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
             values_.resize(nRbf);
             for(int i = 0; i < nRbf; i++){
                 xi = mesh->getXCoord(mesh->getCells(stencil[nStencils-1][i]));
                 xB_[i] = xi;
                 values_[i] = values[nStencils-1][i];
             }
         }
 
         int nRbf0 = std::min(nRbf,3);
         int degPoly_ = std::max(getDegPoly(nRbf),0);
         int nP = multi_index.size();
         int nP_tmp = 0;
         for (int i = 0; i < nP; i++) {
             if (multi_index[i][0] + multi_index[i][1] <= degPoly_) {
                 nP_tmp++;
             }
         }
         nP = nP_tmp;
 
         int m = nRbf + nP;
 
         int approxOrder = std::min(degPoly_ + 1, order);
         getInterpCoeff(xB_, values_, eps, nP, approxOrder);
         values_tmp = values_;
         if (nRbf >1){
             for ( int i = 0; i < nRbf_old; i++){
                 values_tmp[i] -= values_old[i];
             }
 
             for ( int i = 0; i < nP_old; i++){
                 values_tmp[i + nRbf] -= values_old[i + nRbf_old];
             }
         }
         values_old = values_;
         std::vector<std::vector<double>> EvA0(nEv, std::vector<double>(nRbf));
         MVEvaluationMatrix(evPoints, xB_, eps, EvA0, approxOrder);
 
         std::vector<std::vector<double>> PA(nEv, std::vector<double>(nP));
 
         MVPolynomialEvaluationMatrixV(evPoints, xB_, PA, eps, nP);
 
         double *EvA = new double[nEv * m];
 
         for (int i = 0; i < nEv; i++) {
             for (int j = 0; j < nRbf; j++) {
                 EvA[i + j * nEv] = EvA0[i][j];
             }
             for (int j = 0; j < nP; j++) {
                 EvA[i + (j + nRbf) * nEv] = PA[i][j];
             }
         }
 
         if (nRbf > 1){
             smInd[degPoly_] = my_pow(ptwiseSmInd(xB_[0], xB_, values_tmp, eps),1);
         }else{
             smInd[0] = smInd0;
         }
         double *sVecPt = values_tmp.data();
         double *y = new double[nEv];
         cblas_dgemv(CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans, (int) nEv, (int) (nRbf + nP), 1.0, EvA, (int) nEv,
                     sVecPt, 1, 0.0, y, (int) 1);
 
 
         for (int i = 0; i < nEv; i++){
             ret[degPoly_][i] = y[i];
         }
 
         nP_old = nP;
         nRbf_old = nRbf;
 
         delete[] y;
 
         delete[] EvA;
     }
     smInd[0] += 1;
     for (int degPoly_ = 1; degPoly_ <= degPolyMax; degPoly_++){
         smInd[degPoly_] = 1/my_pow(eps0 + my_pow(100,degPoly_)*smInd[degPoly_] + smInd[0],2);
     }
     smInd[0] = 1/my_pow(eps0 + smInd[0] ,2);
 
     for (int degPoly_ = 0; degPoly_ <= degPolyMax; degPoly_++){
         smInd[degPolyMax-degPoly_] /=(smInd[0]);
     }
 
     for ( int j = 0; j <= degPolyMax; j++){
         for (int i = 0; i < nEv; i++){
             res[i] += ret[j][i]*smInd[j];
         }
     }
 
     return res;
 }
 
 std::vector<double> Multiquadratics_2D_WENO_c::evaluateDirect(std::vector<std::vector<double>> evPoints,
                                                             std::vector<std::vector<std::vector<double>>> xB,
                                                             std::vector<double> values, double eps, double rad){
     std::vector<double>res(0);
     return res;
 }
 
 double Multiquadratics_2D_WENO_c::ptwiseSmInd(std::vector<std::vector<double>> xBi,
                                             std::vector<std::vector<std::vector<double>>> xB, std::vector<double> values,
                                             double eps){
 
     int nRbf = (int) xB.size();
     int degPoly_ = std::max(getDegPoly(nRbf),0);
 
     double Ind = 0;
 //    double s = 0;
 //    double p = 1;
 //    double dist = 0;
     std::vector<double> xM(2, 0.0);
 
     xM[0] = (xBi[0][0] + xBi[1][0] + xBi[2][0])/NODE_NUM;
     xM[1] = (xBi[0][1] + xBi[1][1] + xBi[2][1])/NODE_NUM;
 
 //    double delta;
 //    double d_0, d_1;
 //    for (int i = 0; i < nRbf;i++){
 //        Ind += std::pow(values[i],2);
 //        delta = eps*std::sqrt(calcSize(xB[i]));
 //        Ind += values[i]*values[i];
 //        p *= delta*delta;
 //        s += delta*delta;
 //        d_0 = xM[0]-calcMidPt(xB[i],0);
 //        d_1 = xM[1]-calcMidPt(xB[i],1);
 //        dist += eps*std::sqrt(d_0*d_0 + d_1*d_1);
 //    }
 //    Ind = Ind*p*p/(s*s);
 //    Ind = 0;
 
 //    double delta0 = std::sqrt(calcSize(xBi));
     double delta0 = (calcSize(xBi));
     Ind*= delta0*delta0*0;
     for (int j = 1;j <= degPoly_; j++){
 //        double Ds = 0;
 //        for (int i = 0; i < nRbf;i++){
 //            Ds += values[i]*DF(std::vector<double>{eps*(xM[0]-calcMidPt(xB[i],0)),eps*(xM[1]-calcMidPt(xB[i],1))},j);
 //        }
 //        Ds += DPolynomial(xM, xB, values, j, values.size() - nRbf, eps);
 //        Ind += std::pow(Ds,2);
         for (int j2 = 0; j2 <= j; j2++){
             double Ds = 0;
             double Ds_tmp = 0;
             for (int i = 0; i < nRbf;i++){
                 Ds += values[i]*DF(std::vector<double>{(xM[0]-calcMidPt(xB[i],0)),(xM[1]-calcMidPt(xB[i],1))},j2, j-j2, eps);
 //                Ds += values[i]*DF(std::vector<double>{eps*(xM[0]),eps*(xM[1])},j2, j-j2, eps);
 //                Ds_tmp += values[i]*DF(std::vector<double>{eps*(xM[0]-calcMidPt(xB[i],0)),eps*(xM[1]-calcMidPt(xB[i],1))},j2, j-j2, eps);
             }
 //            Ds += DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps);
             double Dt = DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps);
 //            Ind += (std::pow(my_pow(eps,j)*Ds + Dt,2)) *my_pow(delta0,j);
             Ind += (std::pow(my_pow(1,j)*Ds,2) + std::pow(Dt,2)) *my_pow(delta0,j);
         }
     }
 //    if (Ind < 0.1){
 //        std::cout<< Ind;
 //    }
 //    for (int j2 = 0; j2 <= degPoly_; j2++){
 //        double Ds = 0;
 //        for (int i = 0; i < nRbf;i++){
 //                Ds += values[i]*DF(std::vector<double>{eps*(xM[0]),eps*(xM[1])},j2, degPoly_-j2, eps);
 //        }
 ////        Ds += DPolynomial(xM, xB, values, j2, degPoly_-j2, values.size() - nRbf, eps);
 ////            double Dt = DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps);
 //        Ind += std::pow(Ds,2)*my_pow(delta0,degPoly_);
 //    }
 //    Ind *= dist;
 
     return std::fabs(Ind);
 }
 
 int Multiquadratics_2D_WENO_c::getDegPoly(int n){
     int degPoly_;
     if (false){
         degPoly_ = 1;
     }else{
         degPoly_ = (int)floor(-1.5+sqrt(1+8*(n))/2);
     }
     return degPoly_;
 }