diff --git a/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D.cpp b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D.cpp
index ebc28f3..8397ec4 100644
--- a/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D.cpp
+++ b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D.cpp
@@ -1,846 +1,846 @@
 //
 // Created by Fabian Moenkeberg on 31.07.18.
 //
 
 #include "Multiquadratics_2D.hpp"
 #include <numeric>
 
 #include "lapacke.h"
 #include "Utility/triangle_dunavant_rule.hpp"
 #include "Utility/JBurkardt/triangle_wandzura_rule.hpp"
 #include <assert.h>
 #include "../../Utility/vvra.hpp"
 #include "cblas.h"
 #include "Utility/JBurkardt/triangle_integrals.hpp"
 #include <functional>
 #include "Utility/JBurkardt/condition.hpp"
 //#include "Utility/reference_to_physical_t4.cpp"
 
 //#include "Utility/JBurkardt/r8lib.hpp"
 
 Multiquadratics_2D::Multiquadratics_2D(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+3);
 }
 
 Multiquadratics_2D::Multiquadratics_2D(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);
 }
 
 void Multiquadratics_2D::initialize(Configuration config) {
     order_num = dunavant_order_num(order_scheme);
     xy_dunavant = new double[2*order_num];
     weights_dunavant = new double[order_num];
     dunavant_rule(order_scheme, order_num, xy_dunavant, weights_dunavant);
 
     // Quadrilateral weights
     int n = order_scheme;
     std::vector<double> w0(n/2 + (n%2));
     std::vector<double> abs(n/2 + (n%2));
     x_legendre.resize(n);
     weights_legendre.resize(n);
     double *w = new double[n];
     double *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];
         }
     }
 }
 
 void Multiquadratics_2D::setMulti_Index(int n){
     multi_index.resize(n*n,std::vector<int>(2));
     int i0 = 0;
     for (int i = 0; i < n; i++){
         for (int j = 0; j <= i ; j++){
             multi_index[i0][0] = j;
             multi_index[i0][1] = i-j;
             i0++;
         }
     }
 
     int nn = n*n;
 
     for (int j = i0; j < nn; j++){
         multi_index.erase(multi_index.begin()+i0);
     }
 }
 
 
 std::vector<double> Multiquadratics_2D::evaluate(std::vector<std::vector<double>> evPoints, MeshBase *mesh, std::vector<int> stencil, std::vector<double> values){
     int nStencils = (int) stencil.size();
     std::vector<std::vector<std::vector<double>>> xB(nStencils,std::vector<std::vector<double>>(NODE_NUM, std::vector<double>(2)));
     std::vector<std::vector<double>> xi(NODE_NUM,std::vector<double>(2));
     for (int i = 0; i < nStencils; i++){
         xi = mesh->getXCoord(mesh->getCells(stencil[i]));
 //        xB[i].resize(NODE_NUM,std::vector<double>(2));
         xB[i] = xi;
 //        for (int j = 0; j < nNodes; j++){
 //            xB[i][j][0] = xi[j][0];
 //            xB[i][j][1] = xi[j][1];
 //        }
 //        std::cout<< "tri: " << xi[0][0] << "," << xi[0][1] << ","<< xi[1][0] << ","<< xi[1][1] << ","<< xi[2][0] << ","<< xi[2][1] << "\n";
     }
 
     double rad = stableRad / mesh->getApproxStencilSize(stencil);
 //    double rad = stableRad / sqrt(mesh->getSize(cellIndex));
 
     std::vector<double> ret = Multiquadratics_2D::evaluateDirect(evPoints, xB, values, 1, rad);
     return ret;
 };
 
 
 std::vector<double> Multiquadratics_2D::evaluateDirect(std::vector<std::vector<double>> evPoints,
                                                                      std::vector<std::vector<std::vector<double>>> xB,
                                                                      std::vector<double> values, double eps, double rad){
     int nRbf = (int) xB.size();
     int degPoly_ = std::max(getDegPoly(nRbf),0);
 //    degPoly_ = degPoly;
     int nP = multi_index.size();
     int nP_tmp = 0;
     int sum_multi_index = multi_index[0][0]+ multi_index[0][1];
     for (int i = 0; i < nP; i ++){
         if(multi_index[i][0]+ multi_index[i][1] <= degPoly_){
             nP_tmp++;
         }
     }
     nP = nP_tmp;
 //    nP = multi_index.size();
 
     int m = nRbf + nP;
     int nEv = (int) evPoints.size();
     eps = eps*rad;
 
     int approxOrder = std::min(degPoly_+1,order);
     getInterpCoeff(xB, values, eps, nP, approxOrder);
     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];
         }
     }
 
     double* sVecPt = values.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);
 
 
     std::vector<double> ret(y, y+nEv);
     delete[] y;
     delete[] EvA;
 //    delete[] sVecPt;
     return ret;
 }
 
 
 double Multiquadratics_2D::smoothnessInd(int cellIndex, std::vector<int> stencil, std::vector<double> values, MeshBase *mesh){
 
     int nStencils = (int) stencil.size();
     std::vector<std::vector<std::vector<double>>> xB(nStencils,std::vector<std::vector<double>>(NODE_NUM,std::vector<double>(mesh->getNdims())));
     std::vector<std::vector<double>> xi(NODE_NUM,std::vector<double>(2));
 //    int nNodes;
     for (int i = 0; i < nStencils; i++){
         xi = mesh->getXCoord(mesh->getCells(stencil[i]));
 //        nNodes = xi.size();
 //        xB[i].resize(NODE_NUM,std::vector<double>(2));
         xB[i] = xi;
 //        for (int j = 0; j < nNodes; j++){
 //            xB[i][j][0] = xi[j][0];
 //            xB[i][j][1] = xi[j][1];
 //        }
 //        std::cout<< "tri: " << xi[0][0] << "," << xi[0][1] << ","<< xi[1][0] << ","<< xi[1][1] << ","<< xi[2][0] << ","<< xi[2][1] << "\n";
     }
     std::vector<int> evCell = mesh->getCells(cellIndex);
 //    nNodes = evCell.size();
     std::vector<std::vector<double>> xBi(NODE_NUM,std::vector<double>(2));
 
     xBi[0] = mesh->getXCoord(evCell[0]);
     xBi[1] = mesh->getXCoord(evCell[1]);
     xBi[2] = mesh->getXCoord(evCell[2]);
 
     double rad = 1 / sqrt(mesh->getSize(cellIndex));
 
     return ptwiseSmInd(xBi, xB, values, rad);
 //    return 0;
 }
 
 
 void Multiquadratics_2D::MVInterpolationMatrix(std::vector<std::vector<std::vector<double>>> x, double eps, double* M, int order_){
     int nRbf = (int) x.size();
     std::vector<std::vector<double>> node_xy(nRbf, std::vector<double>(2*NODE_NUM));
 
     std::vector<std::vector<double>> xy(nRbf, std::vector<double>(2*order_num));
 //    std::vector<std::vector<double>> xM(nRbf,std::vector<double>(2,0));
     int nNodes;
 //    std::vector<std::vector<double>> M_test(nRbf, std::vector<double>(nRbf));
     for (int i = 0; i < nRbf; i++){
         node_xy[i][0] = x[i][0][0];
         node_xy[i][1] = x[i][0][1];
         node_xy[i][2] = x[i][1][0];
         node_xy[i][3] = x[i][1][1];
         node_xy[i][4] = x[i][2][0];
         node_xy[i][5] = x[i][2][1];
 
         reference_to_physical_t3 ( node_xy[i], order_num, xy_dunavant, xy[i] );
     }
 
     double F_ev;
     double d_0, d_1;
     double eps2 = eps*eps;
 //    eps2 = 0.1;
     for (int i = 0; i < nRbf; i++){
         for (int j = 0; j < nRbf; j++){
             M[j + i*nRbf] = 0;
             for (int k = 0; k<order_num; k++){
                 for (int l = 0; l<order_num; l++) {
                     d_0 = xy[i][0 + k * 2] - xy[j][0 + l * 2];
                     d_1 = xy[i][1 + k * 2] - xy[j][1 + l * 2];
                     F_ev = phi(eps2 * (d_0*d_0 + d_1*d_1), order_);
                     M[j + i*nRbf] += weights_dunavant[k]*weights_dunavant[l] * F_ev;//Note: we do not need to divide by the cell size, since the Jacobian would be the area.
                 }
             }
 //            d_0 = xM[i][0] - xM[j][0];
 //            d_1 = xM[i][1] - xM[j][1];
 //            M_test[i][j] = phi(eps2 * (d_0*d_0 + d_1*d_1), order_);
         }
     }
 }
 
 
 void Multiquadratics_2D::MVInterpolationMatrix(std::vector<std::vector<std::vector<double>>> x, double eps, std::vector<double> &M, int order_){
     int nRbf = (int) x.size();
     std::vector<std::vector<double>> node_xy(nRbf, std::vector<double>(2*NODE_NUM));
 
     std::vector<std::vector<double>> xy(nRbf, std::vector<double>(2*order_num));
 //    std::vector<std::vector<double>> xM(nRbf,std::vector<double>(2,0));
 //    std::vector<std::vector<double>> M_test(nRbf, std::vector<double>(nRbf));
     for (int i = 0; i < nRbf; i++){
         node_xy[i][0] = x[i][0][0];
         node_xy[i][1] = x[i][0][1];
         node_xy[i][2] = x[i][1][0];
         node_xy[i][3] = x[i][1][1];
         node_xy[i][4] = x[i][2][0];
         node_xy[i][5] = x[i][2][1];
 
         reference_to_physical_t3 ( node_xy[i], order_num, xy_dunavant, xy[i] );
     }
 
     double F_ev;
     double d_0, d_1;
     double eps2 = eps*eps;
 //    eps2 = 0.1;
     for (int i = 0; i < nRbf; i++){
         for (int j = 0; j < nRbf; j++){
             M[j + i*nRbf] = 0;
             for (int k = 0; k<order_num; k++){
                 for (int l = 0; l<order_num; l++) {
                     d_0 = xy[i][0 + k * 2] - xy[j][0 + l * 2];
                     d_1 = xy[i][1 + k * 2] - xy[j][1 + l * 2];
                     F_ev = phi(eps2 * (d_0*d_0 + d_1*d_1), order_);
                     M[j + i*nRbf] += weights_dunavant[k]*weights_dunavant[l] * F_ev;//Note: we do not need to divide by the cell size, since the Jacobian would be the area.
                 }
             }
 //            d_0 = xM[i][0] - xM[j][0];
 //            d_1 = xM[i][1] - xM[j][1];
 //            M_test[i][j] = phi(eps2 * (d_0*d_0 + d_1*d_1), order_);
         }
     }
 }
 
 
 void Multiquadratics_2D::MVInterpolationMatrix_ptwise(std::vector<std::vector<std::vector<double>>> x, double eps, std::vector<double> &M, int order_){
     int nRbf = (int) x.size();
 
     std::vector<std::vector<double>> xM(nRbf,std::vector<double>(2,0));
     for (int i = 0; i < nRbf; i++){
         xM[i][0] = (x[i][0][0] + x[i][1][0] + x[i][2][0])/NODE_NUM;
         xM[i][1] = (x[i][0][1] + x[i][1][1] + x[i][2][1])/NODE_NUM;
     }
 
     double d_0, d_1;
     double eps2 = eps*eps;
     for (int i = 0; i < nRbf; i++){
         for (int j = 0; j < nRbf; j++){
             d_0 = xM[i][0] - xM[j][0];
             d_1 = xM[i][1] - xM[j][1];
             M[j + i*nRbf] = phi(eps2 * (d_0*d_0 + d_1*d_1), order_);
         }
     }
 }
 
 
 void Multiquadratics_2D::MVEvaluationMatrix(std::vector<std::vector<double>> x, std::vector<std::vector<std::vector<double>> > xB, double eps, std::vector<std::vector<double>> &EvA, int order_){
 
     int nRbf = (int) xB.size();
     int nEv = (int) x.size();
     std::vector<std::vector<double>> node_xy(nRbf, std::vector<double>(2*NODE_NUM));
 
     std::vector<std::vector<double>> xy(nRbf, std::vector<double>(2*order_num));
     int nNodes;
     for (int i = 0; i < nRbf; i++){
         node_xy[i][0] = xB[i][0][0];
         node_xy[i][1] = xB[i][0][1];
         node_xy[i][2] = xB[i][1][0];
         node_xy[i][3] = xB[i][1][1];
         node_xy[i][4] = xB[i][2][0];
         node_xy[i][5] = xB[i][2][1];
 
         reference_to_physical_t3 ( node_xy[i], order_num, xy_dunavant, xy[i] );
 //        std::cout << xy[i][0] << ", " << xy[i][1]<< ", " <<  xy[i][2] << ", " << xy[i][3]<< ", " << xy[i][4] << ", " << xy[i][5] <<"\n";
     }
 
 //    std::cout << "eps: " << eps << "\n";
 //    std::cout << "xy_dunav: " << xy_dunavant[0] << ", "<< xy_dunavant[1] << ", "<< xy_dunavant[2] << ", "<< xy_dunavant[3] << ", "<< xy_dunavant[4] << ", "<< xy_dunavant[5] << "\n";
 //    for (int i = 0; i < nEv; i++){
 //        std::cout << x[i][0] << ", " << x[i][1] << ",, ";
 //    }
 //    std::cout << "\n";
     double F_ev;
     double d_0;
     double d_1;
     double eps2 = eps*eps;
     for (int i = 0; i < nRbf; i++){
         for (int j = 0; j < nEv; j++){
             EvA[j][i] = 0;
             for (int k = 0; k<order_num; k++){
                 d_0 = xy[i][0+k*2]-x[j][0];
                 d_1 = xy[i][1+k*2]-x[j][1];
                 F_ev = phi(eps2 * (d_0*d_0 + d_1*d_1), order_);
                 EvA[j][i] += weights_dunavant[k]*F_ev; //Note: we do not need to divide by the cell size, since the Jacobian would be the area.
             }
         }
     }
 }
 
 
 void Multiquadratics_2D::MVPolynomialMatrix(std::vector<std::vector<std::vector<double>>> x, double* P, int nm, double eps){
 
 
     int nRbf = (int) x.size();
 
     double *t= new double[2*NODE_NUM];
     for (int i = 0; i < nRbf; i++){
 
         t[0] = eps*(x[i][0][0] - x[0][0][0]);
         t[1] = eps*(x[i][0][1] - x[0][0][1]);
         t[2] = eps*(x[i][1][0] - x[0][0][0]);
         t[3] = eps*(x[i][1][1] - x[0][0][1]);
         t[4] = eps*(x[i][2][0] - x[0][0][0]);
         t[5] = eps*(x[i][2][1] - x[0][0][1]);
 
         for (int j = 0; j <nm; j++) {
             P[i+j*nRbf] = triangle_monomial_average(multi_index[j][0], multi_index[j][1], t);
         }
     }
 
     delete [] t;
 }
 
 
 void Multiquadratics_2D::MVPolynomialMatrix(std::vector<std::vector<std::vector<double>>> x, std::vector<double> &P, int nm, double eps){
 
     int nRbf = (int) x.size();
 
     double *t= new double[2*NODE_NUM];
     for (int i = 0; i < nRbf; i++){
         t[0] = eps*(x[i][0][0] - x[0][0][0]);
         t[1] = eps*(x[i][0][1] - x[0][0][1]);
         t[2] = eps*(x[i][1][0] - x[0][0][0]);
         t[3] = eps*(x[i][1][1] - x[0][0][1]);
         t[4] = eps*(x[i][2][0] - x[0][0][0]);
         t[5] = eps*(x[i][2][1] - x[0][0][1]);
 
         for (int j = 0; j <nm; j++) {
             P[i+j*nRbf] = triangle_monomial_average(multi_index[j][0], multi_index[j][1], t);
         }
     }
 
     delete [] t;
 }
 
 
 void Multiquadratics_2D::MVPolynomialMatrix_ptwise(std::vector<std::vector<std::vector<double>>> x, std::vector<double> &P, int nm, double eps){
 
 
     int nRbf = (int) x.size();
     std::vector<std::vector<double>> xM(nRbf,std::vector<double>(2,0));
     for (int i = 0; i < nRbf; i++){
 
         xM[i][0] = eps*((x[i][0][0] + x[i][1][0] + x[i][2][0])/NODE_NUM - x[0][0][0]);
         xM[i][1] = eps*((x[i][0][1] + x[i][1][1] + x[i][2][1])/NODE_NUM - x[0][0][1]);
 
         for (int j = 0; j <nm; j++) {
             P[i+j*nRbf] = my_pow(xM[i][0], multi_index[j][0])*my_pow(xM[i][1], multi_index[j][1]);
         }
     }
 }
 
 
 void Multiquadratics_2D::MVPolynomialEvaluationMatrixV(std::vector<std::vector<double>> evPoits,
                                                        std::vector<std::vector<std::vector<double>>> xB, std::vector<std::vector<double>> &P, double eps, int nP){
     int nEv = (int)evPoits.size();
 
     for (int i = 0; i < nP; i++){
         for (int j = 0; j < nEv; j++){
             P[j][i] = my_pow(eps*(evPoits[j][0] - xB[0][0][0]),multi_index[i][0])*my_pow(eps*(evPoits[j][1] - xB[0][0][1]),multi_index[i][1]);
         }
     }
 }
 
 
 
 
 void Multiquadratics_2D::getInterpCoeff(std::vector<std::vector<std::vector<double>>> xB, std::vector<double> &values, double eps, int nP, int order_){
     int nRbf = (int) xB.size();
 
     lapack_int info, lda, ldb, nrhs, n, m;
 
     m = nRbf + nP;
     n = m;
     double* s = new double[m];
     int N = (nP+nRbf)*(nP+nRbf);
     double* A = new double[N];
     std::vector<double>A0(nRbf*nRbf);
     MVInterpolationMatrix(xB, eps, A0, order_);
 //    double* A0 = new double[nRbf*nRbf];
 //    MVInterpolationMatrix(xB, eps, A0, order_);
     std::vector<double> PT(nRbf*(nP));
     MVPolynomialMatrix(xB, PT, nP, eps);
 //    double* PT = new double[nRbf*(nP)];
 //    MVPolynomialMatrix(xB, PT, nP, eps);
     for (int i = 0; i < nRbf; i++){
         s[i] = values[i];
         for (int j1 = 0; j1 < nRbf; j1++){
             A[j1 + i*(m)] = A0[j1 + i*(nRbf)];
         }
         for (int j2 = 0; j2 < (nP); j2++){
             A[nRbf + j2 + i*(m)] = PT[i + j2*(nRbf)];
         }
     }
 
     for (int i = 0; i < nP; i++){
         s[i + nRbf] = 0;
         for(int j = 0; j < nRbf; j++){
             A[j + i*m + nRbf*(m)] = PT[j + i*(nRbf)];
         }
         for (int j2 = 0; j2 < nP; j2++){
             A[nRbf + i*(m) + nRbf*(m) + j2] = 0;
         }
     }
 //    std::cout << "\n" << m<<"\n";
 
 //    std::vector<double> testA(A, A + N);
 //    double cond = condition_hager ( nP+nRbf, A );
 //
 //    if (cond > 10000000)
 //        std::cout << cond << "\n";
     nrhs = 1;
     lda = n;
     ldb = m;
     info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m,n,nrhs,A,lda,s,ldb);
 
     values.insert(values.begin(),s,s+m);
     values.resize(m);
 
 //    delete [] PT;
 //    delete [] A0;
     delete [] A;
     delete [] s;
 }
 
 
 void Multiquadratics_2D::getInterpCoeff_ptwise(std::vector<std::vector<std::vector<double>>> xB, std::vector<double> &values, double eps, int nP, int order_){
     int nRbf = (int) xB.size();
 
     lapack_int info, lda, ldb, nrhs, n, m;
 
     m = nRbf + nP;
 //    std::cout << nRbf << ", " << nP<< ", " << values.size() << "\n";
     n = m;
     double* s = new double[m];
     int N = (nP+nRbf)*(nP+nRbf);
     double* A = new double[N];
     std::vector<double>A0(nRbf*nRbf);
     MVInterpolationMatrix_ptwise(xB, eps, A0, order_);
 
     std::vector<double> PT(nRbf*(nP));
     MVPolynomialMatrix_ptwise(xB, PT, nP, eps);
     for (int i = 0; i < nRbf; i++){
         s[i] = values[i];
         for (int j1 = 0; j1 < nRbf; j1++){
             A[j1 + i*(m)] = A0[j1 + i*(nRbf)];
         }
         for (int j2 = 0; j2 < (nP); j2++){
             A[nRbf + j2 + i*(m)] = PT[i + j2*(nRbf)];
         }
     }
 
     for (int i = 0; i < nP; i++){
         s[i + nRbf] = 0;
         for(int j = 0; j < nRbf; j++){
             A[j + i*m + nRbf*(m)] = PT[j + i*(nRbf)];
         }
         for (int j2 = 0; j2 < nP; j2++){
             A[nRbf + i*(m) + nRbf*(m) + j2] = 0;
         }
     }
 //    std::cout << "\n" << m<<"\n";
 
 //    std::vector<double> testA(A, A + N);
 //    double cond = condition_hager ( nP+nRbf, A );
 //
 //    if (cond > 10000000)
 //        std::cout << cond << "\n";
     nrhs = 1;
     lda = n;
     ldb = m;
     info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m,n,nrhs,A,lda,s,ldb);
 
     values.insert(values.begin(),s,s+m);
     values.resize(m);
 
 //    delete [] PT;
 //    delete [] A0;
     delete [] A;
     delete [] s;
 }
 
 
 double Multiquadratics_2D::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);
 
     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;
 
 //    getInterpCoeff(xB, values, eps, nP, order);
     int approxOrder = std::min(degPoly_+1,order);
 //    getInterpCoeff(xB, values, eps, nP, approxOrder);
     getInterpCoeff_ptwise(xB, values, eps, nP, approxOrder);
     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 = eps*std::sqrt(calcSize(xBi));
     for (int j = 1;j <= degPoly_+1; 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;
 //            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))},j2, j-j2, eps);
 //            }
             Ds += DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps);
             Ind += std::pow(Ds,2)*my_pow(delta0,j);
         }
     }
 //    Ind *= dist;
 
     return std::fabs(Ind);
 }
 
 
 double Multiquadratics_2D::phi(double x2) {
 //    return std::pow((1.0 + x2),(order-0.5));
 //    return std::sqrt(1.0 + x2);
     x2 = 1.0 + x2;
     return my_pow(x2, order - 1)*std::sqrt(x2);
 }
 
 double Multiquadratics_2D::phi(double x2, int order_) {
 //    return std::pow((1.0 + x2),(order_-0.5));
 //    return std::sqrt(1.0 + x2);
     x2 = 1.0 + x2;
     return my_pow(x2, order_ - 1)*std::sqrt(x2);
 }
 
 
 double Multiquadratics_2D::DfPhi(double x, int n){
     long double ret;
     switch (n) {
         case 0:
             x = 1.0 + x;
             ret = my_pow(x, order - 1)*std::sqrt(x);
             break;
         case 1:
-            ret = (2*order - 1)/2*pow(1 + x,(2*(double)order - 3)/2);
+            ret = (2*order - 1)/2.0*pow(1 + x,(2*(double)order - 3)/2.0);
 //            x = 1.0 + x;
 //            ret = my_pow(x,order - 2)*std::sqrt(x);
             break;
         case 2:
-            ret = (2*order - 1)*(2*order - 3)/4*pow(1 + x,(2*(double)order - 5)/2);
+            ret = (2*order - 1)*(2*order - 3)/4.0*pow(1 + x,(2*(double)order - 5)/2.0);
 //            x = 1.0 + x;
 //            ret = my_pow(x,order - 3)*std::sqrt(x);
             break;
         case 3:
-            ret = (2*order - 1)*(2*order - 3)*(2*order - 5)/8*pow(1 + x,(2*(double)order - 7)/2);
+            ret = (2*order - 1)*(2*order - 3)*(2*order - 5)/8.0*pow(1 + x,(2*(double)order - 7)/2.0);
 //            x = 1.0 + x;
 //            ret = my_pow(x,order - 4)*std::sqrt(x);
             break;
         case 4:
-            ret = (2*order - 1)*(2*order - 3)*(2*order - 5)*(2*order - 7)/16*pow(1 + x,(2*(double)order - 9)/2);
+            ret = (2*order - 1)*(2*order - 3)*(2*order - 5)*(2*order - 7)/16.0*pow(1 + x,(2*(double)order - 9)/2.0);
 //            x = 1.0 + x;
 //            ret = my_pow(x,order - 5)*std::sqrt(x);
             break;
         case 5:
-            ret = (2*order - 1)*(2*order - 3)*(2*order - 5)*(2*order - 7)*(2*order - 9)/32*pow(1 + x,(2*(double)order - 11)/2);
+            ret = (2*order - 1)*(2*order - 3)*(2*order - 5)*(2*order - 7)*(2*order - 9)/32.0*pow(1 + x,(2*(double)order - 11)/2.0);
 //            x = 1.0 + x;
 //            ret = my_pow(x,order - 6)*std::sqrt(x);
             break;
         default:
             assert(0 && "Multiquadratics (DfPhi): higher orders not implemented.");
             break;
     }
     return ret;
 }
 
 double Multiquadratics_2D::DPolynomial(std::vector<double> x, std::vector<std::vector<std::vector<double>>> xB, std::vector<double> coeff, int i, int nP, double eps){
 //    int nP = multi_index.size();
     int n = coeff.size();
     double ret = 0.0;
 
     for (int j = 0; j < nP; j ++){
         double temp = 0.0;
 //        std::cout << "MultiInd: " << multi_index[j][0] << ", " << multi_index[j][1] << "\n";
         for ( int l = 0; l <= i; l++){
             int temp0 = 1;
             int temp1 = 1;
             double temp2 = 1;
             if (multi_index[j][0]-i+l < 0||multi_index[j][1]-l < 0){
             }
             else {
                 for (int k = 0; k < i - l; k++) {
                     temp0 *= (multi_index[j][0] - k);
 //                std::cout << (multi_index[j][0]-k) << ", " ;
                 }
 //            std::cout << "\n" ;
                 for (int k = 0; k < l; k++) {
                     temp1 *= (multi_index[j][1] - k);
                 }
                 if (l == 0) {
                     temp2 = 1;
                 } else {
                     for (int k = 0; k < i - l; k++) {
                         temp2 *= double((i - k)) / double((l - k));
 //                    std::cout << i - k << ", " << l - k << "\n ";
                     }
                 }
 //                if (temp2>1)
 //                    std::cout << "temp: ";
 //                std::cout << "temp: " << temp1 << ", " << temp2 << "," << temp0 << "\n";
 //                std::cout << "Index: " << multi_index[j][0] - i + l << ", " << multi_index[j][1] - l << "\n";
 //                std::cout << "Index: " << pow(x[0], multi_index[j][0] - i + l) << ", " << pow(x[1], multi_index[j][1] - l) << "," << temp+temp2 * temp0 * temp1 * pow(x[0], multi_index[j][0] - i + l) * pow(x[1], multi_index[j][1] - l) << "\n";
 //                temp += temp2 * temp0 * temp1 * pow(eps*(x[0] - xB[0][0][0]), multi_index[j][0] - i + l) * pow(eps*(x[1] - xB[0][0][1]), multi_index[j][1] - l);
                 temp += temp0 * temp1 * pow(eps*(x[0] - xB[0][0][0]), multi_index[j][0] - i + l) * pow(eps*(x[1] - xB[0][0][1]), multi_index[j][1] - l);
             }
         }
 //        temp *= pow(eps,i);
         ret += coeff[n-nP+j]*temp;
     }
     return ret;
 }
 
 double Multiquadratics_2D::DPolynomial(std::vector<double> x, std::vector<std::vector<std::vector<double>>> xB, std::vector<double> coeff, int i, int i2, int nP, double eps){
 //    int nP = multi_index.size();
     int n = coeff.size();
     double ret = 0.0;
 
     for (int j = 0; j < nP; j ++){
         double temp = 0.0;
         int temp0 = 1;
         int temp1 = 1;
         if (multi_index[j][0]-i2 < 0||multi_index[j][1]-i < 0){
         }
         else {
             for (int k = 0; k < i2; k++) {
                 temp0 *= (multi_index[j][0] - k);
             }
             for (int k = 0; k < i; k++) {
                 temp1 *= (multi_index[j][1] - k);
             }
 
             temp += temp0 * temp1 * pow(eps*(x[0] - xB[0][0][0]), multi_index[j][0] - i2) * pow(eps*(x[1] - xB[0][0][1]), multi_index[j][1] - i);
         }
 //        temp *= pow(eps,i);
         ret += coeff[n-nP+j]*temp;
     }
     return ret;
 }
 
 
 double Multiquadratics_2D::DF(std::vector<double> x, int i){
     long double ret;
     switch (i) {
         case 0:
             ret = DfPhi(x[0]*x[0]+x[1]*x[1],0);
             break;
         case 1:
             ret = 2*(x[0]+x[1])*DfPhi(x[0]*x[0]+x[1]*x[1],1);
             break;
         case 2:
             ret = 4*(DfPhi(x[0]*x[0]+x[1]*x[1],1)+(x[0]*x[0]+x[0]*x[1]+x[1]*x[1])*DfPhi(x[0]*x[0]+x[1]*x[1],2));
             break;
         case 3:
             ret = 8*(x[0]+x[1])*(2*DfPhi(x[0]*x[0]+x[1]*x[1],2)+(x[0]*x[0]+x[1]*x[1])*DfPhi(x[0]*x[0]+x[1]*x[1],3));
             break;
         case 4:
             ret = 4*(7*DfPhi(x[0]*x[0]+x[1]*x[1],2)+2*(7*x[0]*x[0]+6*x[0]*x[1]+7*x[1]*x[1])*DfPhi(x[0]*x[0]+x[1]*x[1],3)
                     +4*(pow(x[0],4)+pow(x[0],3)*x[1]+x[0]*x[0]*x[1]*x[1]+pow(x[1],3)*x[0]+pow(x[1],4))*DfPhi(x[0]*x[0]+x[1]*x[1],4));
             break;
         case 5:
             ret = 8*(x[0]+x[1])*(21*DfPhi(x[0]*x[0]+x[1]*x[1],3)+(22*x[0]*x[0]-4*x[0]*x[1]+22*x[1]*x[1])*DfPhi(x[0]*x[0]+x[1]*x[1],4)
                      +4*(pow(x[0],4)+x[0]*x[0]*x[1]*x[1]+pow(x[1],4))*DfPhi(x[0]*x[0]+x[1]*x[1],5));
             break;
         default:
             assert(0 && "Multiquadratics (DF): higher orders not implemented.");
             break;
     }
     return ret;
 }
 
 
 double Multiquadratics_2D::DF(std::vector<double> x, int i, int j, double eps){
     long double ret;
     if (i==0 && j==0){
         ret = DfPhi(x[0]*x[0]+x[1]*x[1],0);
     }
     else if (i==1 && j==0){
         ret = 2*x[0]*DfPhi(x[0]*x[0]+x[1]*x[1],1);
     }else if (i==0 && j==1){
         ret = 2*x[1]*DfPhi(x[0]*x[0]+x[1]*x[1],1);
     }else if (i==2 && j==0){
         ret = 2*eps*eps*DfPhi(x[0]*x[0]+x[1]*x[1],1) + 4*x[0]*x[0]*DfPhi(x[0]*x[0]+x[1]*x[1],2);
     }else if (i==1 && j==1){
         ret = 4*x[0]*x[1]*DfPhi(x[0]*x[0]+x[1]*x[1],2);
     }else if (i==0 && j==2){
         ret = 2*eps*eps*DfPhi(x[0]*x[0]+x[1]*x[1],1) + 4*x[1]*x[1]*DfPhi(x[0]*x[0]+x[1]*x[1],2);
     }else if (i==3 && j==0) {
         ret = 12 * eps * x[0] * DfPhi(x[0] * x[0] + x[1] * x[1], 2) + 8 * x[0] * x[0] * x[0] * DfPhi(x[0] * x[0] + x[1] * x[1], 3);
     }else if (i==2 && j==1) {
         ret = 4 * eps * x[1] * DfPhi(x[0] * x[0] + x[1] * x[1], 2) + 8 * x[0] * x[0] * x[1] * DfPhi(x[0] * x[0] + x[1] * x[1], 3);
     }else if (i==1 && j==2) {
         ret = 4 * x[0] * eps * DfPhi(x[0] * x[0] + x[1] * x[1], 2) + 8 * x[1] * x[1] * x[0] * DfPhi(x[0] * x[0] + x[1] * x[1], 3);
     }else if (i==0 && j==3) {
         ret = 12 * x[1] * eps * DfPhi(x[0] * x[0] + x[1] * x[1], 2) + 8 * x[1] * x[1] * x[1] * DfPhi(x[0] * x[0] + x[1] * x[1], 3);
     }else{
         assert(0 && "Multiquadratics (DF): higher orders not implemented.");
     }
     return ret;
 }
 
 
 double Multiquadratics_2D::DF2(std::vector<double> x, int i){
     long double ret;
     switch (i) {
         case 0:
             ret = DfPhi(x[0]*x[0]+x[1]*x[1],0);
             break;
         case 1:
             ret = 4*(x[0]*x[0]+x[1]*x[1])*pow(DfPhi(x[0]*x[0]+x[1]*x[1],1),2);
             break;
         case 2:
             ret = 8*(pow(DfPhi(x[0]*x[0]+x[1]*x[1],1),2)+2*(x[0]*x[0]+x[1]*x[1])*DfPhi(x[0]*x[0]+x[1]*x[1],1)*DfPhi(x[0]*x[0]+x[1]*x[1],2)+2*(pow(x[0],4)+pow(x[0],2)*pow(x[1],2)+pow(x[1],4))*pow(DfPhi(x[0]*x[0]+x[1]*x[1],2),2));
             break;
         case 3:
             ret = 32*(5*(x[0]*x[0]+x[1]*x[1])*pow(DfPhi(x[0]*x[0]+x[1]*x[1],1),2)+2*(3*pow(x[0],4)+2*x[0]*x[0]*x[1]*x[1]+3*pow(x[1],4))*DfPhi(x[0]*x[0]+x[1]*x[1],2)*DfPhi(x[0]*x[0]+x[1]*x[1],3)
                     +2*(pow(x[0],6)+pow(x[0],4)*pow(x[1],2)+pow(x[0],2)*pow(x[1],4)+pow(x[1],6))*pow(DfPhi(x[0]*x[0]+x[1]*x[1],3),2));
             break;
         default:
             assert(0 && "Multiquadratics (DF): higher orders not implemented.");
             break;
     }
     return ret;
 }
 
 double Multiquadratics_2D::my_pow(double x, int n){
     double r = 1.0;
 
     while(n > 0){
         r *= x;
         --n;
     }
 
     return r;
 }
 
 int Multiquadratics_2D::getDegPoly(int n){
     int degPoly_;
     if (n>9){
         degPoly_ = 2;
     }else{
 //        degPoly_ = floor(-1.5+sqrt(1+8*n)/2);
 //        degPoly_ = (int)floor(-1.5+sqrt(1+8*(n-1))/2);
         degPoly_ = (int)floor(-1.5+sqrt(1+8*(n))/2);
     }
     return degPoly_;
 }
 
 double Multiquadratics_2D::calcSize(std::vector<std::vector<double>> xBi){
     return 0.5*std::fabs((xBi[0][0]-xBi[2][0])*(xBi[1][1]-xBi[0][1])-(xBi[0][0]-xBi[1][0])*(xBi[2][1]-xBi[0][1]));
 }
 
 double Multiquadratics_2D::calcMidPt(std::vector<std::vector<double>> xBi, int i){
 
     return (xBi[0][i] + xBi[0][1] + xBi[2][i])/NODE_NUM;
 }
\ No newline at end of file