diff --git a/ConSol/ConSol/Base/Configuration.cpp b/ConSol/ConSol/Base/Configuration.cpp
index 0dff95e..24fb6ef 100644
--- a/ConSol/ConSol/Base/Configuration.cpp
+++ b/ConSol/ConSol/Base/Configuration.cpp
@@ -1,883 +1,885 @@
 //
 //  Configuration.cpp
 //  ConserSolver
 //
 //  Created by Fabian Moenkeberg on 04.05.17.
 //  Copyright © 2017 FabianMoenkeberg. All rights reserved.
 //
 
 #include "Configuration.hpp"
 #include <omp.h>
 #include <iomanip>
 
 Configuration::Configuration(){
 };
 
 Configuration::Configuration(char* file) {
 
     this->input_file = file;
     read();
     this->SOLN_DIR = OUTPUT_PATH + "/SOLN_FILES";
 
 //    // Making solution directories
 //    std::string commandline0 = "rm -rf " + SOLN_DIR;
 //    system(commandline0.c_str());
 //    std::string commandline1 = "mkdir -p " + SOLN_DIR;
 //    system(commandline1.c_str());
 //
 //    glob_filename = SOLN_DIR + "/globals.dat";
 //    glob_file.open (glob_filename.c_str());
 //
 //    L1_err_filename = SOLN_DIR + "/l1_err.dat";
 //    L2_err_filename = SOLN_DIR + "/l2_err.dat";
 //    Linf_err_filename = SOLN_DIR + "/linf_err.dat";
 //
 //    if(param.exact_available)
 //    {
 //        L1_err_file.open (L1_err_filename.c_str());
 //        L2_err_file.open (L2_err_filename.c_str());
 //        Linf_err_file.open (Linf_err_filename.c_str());
 //    }
 //
 //    if(param.OOC)
 //    {
 //        print_globals = false;
 //        print_soln = false;
 //        disp_log = false;
 //        L1_err.resize(param.mesh_levels);
 //        L2_err.resize(param.mesh_levels);
 //        Linf_err.resize(param.mesh_levels);
 //    }
 //    else {
 //        print_globals = true;
 //        print_soln = true;
 //        disp_log = true;
 //    }
 }
 
 Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, int order, double tMax, double cfl){
     this->model = model;
     this->evalMethod = evalMethod;
     this->bcname = bc;
     this->order = order;
     this->tMax = tMax;
     this->cfl = cfl;
     this->tInclude = *new std::vector<double>(0);
     this->reconstr = NORECONSTR;
     this->interpBase = POLYNOMIALS;
 }
 
 Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, int order, double tMax, std::vector<double> tInclude, double cfl){
     this->model = model;
     this->evalMethod = evalMethod;
     this->bcname = bc;
     this->order = order;
     this->tMax = tMax;
     this->cfl = cfl;
     this->tInclude = tInclude;
     this->reconstr = NORECONSTR;
     this->interpBase = POLYNOMIALS;
 }
 
 Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, int order, double tMax, double cfl){
     this->model = model;
     this->evalMethod = evalMethod;
     this->bcname = bc;
     this->order = order;
     this->tMax = tMax;
     this->cfl = cfl;
     this->tInclude = *new std::vector<double>(0);
     this->reconstr = reconstr;
     this->interpBase = interpBase;
 }
 
 Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, int order, double tMax, std::vector<double> tInclude, double cfl){
     this->model = model;
     this->evalMethod = evalMethod;
     this->bcname = bc;
     this->order = order;
     this->tMax = tMax;
     this->cfl = cfl;
     this->tInclude = tInclude;
     this->reconstr = reconstr;
     this->interpBase = interpBase;
 }
 
 Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, bool reconstrOnCharact, int order, double tMax, double cfl){
     this->model = model;
     this->evalMethod = evalMethod;
     this->bcname = bc;
     this->order = order;
     this->tMax = tMax;
     this->cfl = cfl;
     this->tInclude = *new std::vector<double>(0);
     this->reconstr = reconstr;
     this->interpBase = interpBase;
     this->reconstrOnCharact = reconstrOnCharact;
 }
 
 Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, bool reconstrOnCharact, int order, double tMax, std::vector<double> tInclude, double cfl){
     this->model = model;
     this->evalMethod = evalMethod;
     this->bcname = bc;
     this->order = order;
     this->tMax = tMax;
     this->cfl = cfl;
     this->tInclude = tInclude;
     this->reconstr = reconstr;
     this->interpBase = interpBase;
     this->reconstrOnCharact = reconstrOnCharact;
 }
 
 ModelName Configuration::getModel(){return model;};
 
 TimeIntegrationName Configuration::getTimeIntegration(){return time_scheme;};
 
 int Configuration::getOrder(){return order;};
 
 evaluationMethod Configuration::getEvalMethod(){return evalMethod;};
 
 bcName Configuration::getBc(){return bcname;};
 
 std::vector<double> Configuration::getTInclude(){return tInclude;};
 
 limiterName  Configuration::getLimiter(){return limiter_name; };
 
 FluxName  Configuration::getFlux(){return fluxScheme; };
 
 double Configuration::getTmax(){return tMax;};
 
 double Configuration::getCfl(){return cfl;};
 
 int Configuration::getMaxNumWrite(){return maxNumWrite;};
 
 ReconstrName Configuration::getReconstr(){return reconstr;};
 
 sourceName Configuration::getSource() { return source;};
 
 initalCondName  Configuration::getInitialCond() {return initialCond;};
 
 InterpBaseName Configuration::getInterpBase(){return interpBase;};
 
 bool Configuration::getReconstrOnCharact(){return reconstrOnCharact;};
 
 int Configuration::getDim(){return dim;};
 
 double Configuration::getXmax(){return xmax;};
 
 double Configuration::getXmin(){return xmin;};
 
 double Configuration::getYmax(){return ymax;};
 
 double Configuration::getYmin(){return ymin;};
 
 double Configuration::getShapeParam(){ return shapeParam;};
 
 std::vector<int>  Configuration::getNx() { return Nx_cell;};
 
 std::vector<int>  Configuration::getNy() { return Ny_cell;};
 
 std::string Configuration::getGridFile() { return grid_file;};
 
 std::string Configuration::getOutputPath() { return OUTPUT_PATH;}
 
 std::string Configuration::getImportPath() { return IMPORT_PATH;}
 
 bool Configuration::getImportSol() {return import_sol;};
 
 bool Configuration::getLoadGrid() { return load_grid;};
 
 int Configuration::getNthreads() {return nThreads;};
 
 stencilChoice Configuration::getStencilChoice() { return  stencilchoice;};
 
 std::map<std::string, double> Configuration::getConstants() {return constants;};
 
 int Configuration::getDegPoly() {return degPoly;};
 
 MeshType Configuration::getMeshType() {return meshType;};
 
 //------------------------------------------------------------------------------
 // Read parameters from file
 //------------------------------------------------------------------------------
 void Configuration::read ()
 {
     std::cout << "Reading input file " << input_file << std::endl;
     Reader fin(input_file);
 
     read_grid (fin);
     read_numeric (fin);
     read_model(fin);
     read_constants (fin);
     read_initial_condition (fin);
     read_exact (fin);
     read_boundary_condition (fin);
     read_output (fin);
     read_import(fin);
 }
 
 //------------------------------------------------------------------------------
 // Read grid section
 //------------------------------------------------------------------------------
 void Configuration::read_grid (Reader &fin)
 {
     std::cout << "  Reading grid section\n";
 
     std::string input;
     int N;
 
     fin.begin_section ("grid");
 
     fin.entry ("load_grid");
     fin >> input;
     if (input == "true") {
         load_grid = true;
     }
     else{
         load_grid = false;
     }
 
     fin.entry ("dim");
     fin >> dim;
 
     assert(dim == 1 || dim == 2 );
 
     fin.entry ("x_min");
     fin >> xmin;
     fin.entry ("x_max");
     fin >> xmax;
 
 
     fin.begin_section ("Nx_cell");
     while(!fin.eos())
     {
         fin >> N;
         assert(N > 0);
         Nx_cell.push_back(N);
     }
     if (dim == 2) {
         fin.entry("y_min");
         fin >> ymin;
         fin.entry("y_max");
         fin >> ymax;
     }else{
         ymin = 0;
         ymax = 1;
     }
 
     fin.begin_section ("Ny_cell");
     while(!fin.eos())
     {
         fin >> N;
         assert(N > 0);
         Ny_cell.push_back(N);
     }
     assert(xmin<xmax);
     assert(ymin<ymax);
 
     fin.entry ("grid_path");
     fin >> grid_file;
 
     fin.entry ("type");
     fin >> input;
     if (input == "triangular") {
         meshType = TRI;
     } else if (input == "hybrid") {
         meshType = HYBRID;
     }else if (input == "quadrilateral"){
         meshType = QUAD;
     }else if (input == "mixed"){
         meshType = MIXED;
     } else {
         std::cout << "   Error: unknown meshType " << input << std::endl;
         exit(0);
     }
 
     fin.end_section ();
 }
 
 //------------------------------------------------------------------------------
 // Read numeric section
 //------------------------------------------------------------------------------
 void Configuration::read_numeric (Reader &fin) {
     std::cout << "  Reading numeric section\n";
 
     std::string input;
 
     fin.begin_section("numeric");
 
     fin.entry("nthreads");
     fin >> nThreads;
 
     fin.entry("time_scheme");
     fin >> input;
     if (input == "SSPRK3") {
         time_scheme = SSPRK3;
     } else if (input == "SSPRK5") {
         time_scheme = SSPRK5;
     }else if (input == "SSPRK2"){
         time_scheme = SSPRK2;
     }else if (input == "FE") {
         time_scheme = FE;
     } else {
         std::cout << "   Error: unknown time_scheme " << input << std::endl;
         exit(0);
     }
 
     fin.entry("cfl");
     fin >> cfl;
     assert (cfl > 0.0);
 
     fin.entry("max_iter");
     fin >> max_iter;
     assert (max_iter > 0);
 
     fin.entry("final_time");
     fin >> tMax;
     assert (tMax >= 0.0);
 
     fin.entry("min_residue");
     fin >> min_residue;
     assert (min_residue > 0.0);
 
     fin.entry("reconstruct");
     fin >> input;
     if (input == "noReconstruct") {
         reconstr = NORECONSTR;
     } else if (input == "eno") {
         reconstr = ENO;
     } else if (input == "weno") {
         reconstr = WENO;
     } else {
         std::cout << "   Error: unknown reconstruction method " << input << std::endl;
         exit(0);
     }
 
     fin.entry("stencilChoice");
     fin >> input;
     if (input == "std") {
         stencilchoice = STD;
     } else if (input == "otf") {
         stencilchoice = OTF;
     } else if (input == "symm") {
         stencilchoice = SYMM;
+    } else if (input == "symmb") {
+        stencilchoice = SYMMb;
     } else {
         std::cout << "   Error: unknown stencil choice " << input << std::endl;
         exit(0);
     }
 
     fin.entry("order");
     fin >> order;
 
     fin.entry("degPoly");
     fin >> input;
 
     if (input == "x"){
         degPoly = order - 2;
     } else {
         degPoly = std::stoi(input);
     }
 
     fin.entry("interpolationBase");
     fin >> input;
     if (input == "Multiquadratics") {
         interpBase = MULTIQUADS;
     } else if (input == "Polynomials") {
         interpBase = POLYNOMIALS;
     } else if (input == "Polyharmonics") {
         interpBase = POLYHARMONICS;
     } else if (input == "Gaussians") {
         interpBase = GAUSSIANS;
     } else {
         std::cout << "   Error: unknown interpolation base " << input << std::endl;
         exit(0);
     }
 
     fin.entry("shapeParameter");
     fin >> shapeParam;
 
     fin.entry("reconstrOnCharact");
     fin >> input;
     if (input == "true") {
         reconstrOnCharact = true;
     }
     else{
         reconstrOnCharact = false;
     }
 
     fin.entry("evaluationMethod");
     fin >> input;
     if (input == "meanvalue") {
         evalMethod = MEANVALUE;
     }
     else{
         evalMethod = PTWISE;
     }
 
     fin.entry("fluxScheme");
     fin >> input;
     if (input == "LaxFried") {
         fluxScheme = LAXFRIEDRICHS;
     }
     else {
         std::cout << "   Error: unknown flux scheme " << input << std::endl;
         exit(0);
     }
 
     fin.entry("limiter");
     fin >> input;
     if (input == "noLimiter") {
         limiter_name = NOLIMITER;
     }
     else if (input == "extremaPreserving") {
         limiter_name = EXTREMAPRES;
     }
     else {
         std::cout << "   Error: unknown limiter " << input << std::endl;
         exit(0);
     }
 
     fin.begin_section("tInclude");
     double value;
 
     this->tInclude = std::vector<double>(0);
     while (!fin.eos())
     {
         fin >> value;
         tInclude.push_back(value);
     }
     fin.end_section ();
 }
 
 //------------------------------------------------------------------------------
 // Read material section
 //------------------------------------------------------------------------------
 void Configuration::read_model (Reader &fin)
 {
     std::cout << "  Reading model section\n";
 
     std::string input;
 
     fin.begin_section ("model");
 
     fin.entry("model");
     fin >> input;
     if (input == "linAdv") {
         model = LINEARADVECTION;
     } else if (input == "burgers") {
         model = BURGERS;
     } else if (input == "SW") {
         model = SW;
     } else if (input == "Euler") {
         model = EULER;
     } else if (input == "BuckleyLeverett") {
         model = BUCKLEYLEVERETT;
     }else if (input == "KPP"){
         model = KPP;
     }else if (input == "KPPY"){
         model = KPPY;
     }else if (input == "KPPX"){
         model = KPPX;
     }
     else {
         std::cout << "   Error: unknown model " << input << std::endl;
         exit(0);
     }
 
     fin.entry("source");
     fin >> input;
     if (input == "noSource") {
         source = NOSOURCE;
     }else{
         std::cout << "   Error: unknown source " << input << std::endl;
         exit(0);
     }
     fin.end_section ();
 }
 
 //------------------------------------------------------------------------------
 // Read constants section
 //------------------------------------------------------------------------------
 void Configuration::read_constants (Reader &fin)
 {
     std::cout << "  Reading constants section\n";
 
     std::string input;
     double value;
 
     fin.begin_section ("constants");
 
     while (!fin.eos())
     {
         fin >> input;
         fin >> value;
         std::cout << std::setw(16) << input << std::setw(16) << value << std::endl;
         constants.insert ( std::pair<std::string,double>(input, value) );
     }
 
 }
 
 //------------------------------------------------------------------------------
 // Read initial condition
 //------------------------------------------------------------------------------
 void Configuration::read_initial_condition (Reader &fin)
 {
     std::cout << "  Reading initial condition section\n";
 
     std::string input;
 
     fin.begin_section ("initial_condition");
 
     switch (model){
         case LINEARADVECTION:
             fin.entry("y");
             fin >> input;
 
             if (input == "sin(2*pi*x)"){
                 initialCond = SIN;
             }else if (input == "bubble") {
                 initialCond = BUBBLE;
             }
             else if ( input == "coscos") {
                 initialCond = COSCOS;
             }
             else if ( input == "sinsin") {
                 initialCond = SINSIN;
             }
             else if (input == "shockX"){
                 initialCond = SHOCKX;
             }
             else if (input == "cube"){
                 initialCond = InitCUBE;
             }else {
                 std::cout << "   Error:: unknown initial Condition: " << input << std::endl;
                 exit (0);
             }
 
             break;
 
         case EULER:
             fin.entry("y");
             fin >> input;
 
             if (input == "isentropicVortex"){
                 initialCond = InitISENTROPICVORTEX;
             }else if (input == "eulerBlast"){
                 initialCond = InitEULERBLAST;
             }else if (input == "vortexShock"){
                 initialCond = InitVORTEXSHOCK;
             }else if (input == "RP"){
                 initialCond = RP;
             }else if (input == "initDoubleMachRef"){
                 initialCond = InitDOUBLEMACHREF;
             }else if (input == "initFAC"){
                 initialCond = InitFAC;
             }else if (input == "initAirFoil"){
                 initialCond = InitAIRFOIL;
             }else if (input == "initNozzle"){
                 initialCond = InitNOZZLE;
             }
             else {
                 std::cout << "   Error:: unknown initial Condition: " << input << std::endl;
                 exit (0);
             }
             break;
         case BURGERS:
             fin.entry("y");
             fin >> input;
 
             if (input == "sin(2*pi*x)"){
                 initialCond = SIN;
             }else if (input == "bubble") {
                 initialCond = BUBBLE;}
             else if (input == "burgers") {
                 initialCond = InitBURGERS;}
             else{
                 std::cout << "   Error:: unknown initial Condition: " << input << std::endl;
                 exit (0);
             }
             break;
         case SW:
             std::cout << "   Error:: Not implemented, yet: " << input << std::endl;
             exit (0);
             break;
         case BUCKLEYLEVERETT:
             std::cout << "   Error:: Not implemented, yet: " << input << std::endl;
             exit (0);
             break;
         case KPP:
             fin.entry("y");
             fin >> input;
 
             if (input == "KPP") {
                 initialCond = InitKPP;
             }else {
                 std::cout << "   Error:: unknown initial Condition: " << input << std::endl;
                 exit (0);
             }
             break;
         case KPPY:
             fin.entry("y");
             fin >> input;
 
             if (input == "KPP") {
                 initialCond = InitKPP;
             }else if (input == "KPPY") {
                 initialCond = InitKPPY;
             }else {
                 std::cout << "   Error:: unknown initial Condition: " << input << std::endl;
                 exit (0);
             }
             break;
         case KPPX:
             fin.entry("y");
             fin >> input;
 
             if (input == "KPP") {
                 initialCond = InitKPP;
             }else {
                 std::cout << "   Error:: unknown initial Condition: " << input << std::endl;
                 exit (0);
             }
             break;
     }
 
 //
 //
 //    fin.entry("density");
 //    fin.getline (input);
 //    initial_condition.add ("density", input);
 //
 //    fin.entry ("xvelocity");
 //    fin.getline (input);
 //    initial_condition.add ("xvelocity", input);
 //
 //    fin.entry ("yvelocity");
 //    fin.getline (input);
 //    initial_condition.add ("yvelocity", input);
 //
 //    fin.entry ("zvelocity");
 //    fin.getline (input);
 //    initial_condition.add ("zvelocity", input);
 //
 //    fin.entry ("pressure");
 //    fin.getline (input);
 //    initial_condition.add ("pressure", input);
 //
 //
     fin.end_section ();
 }
 
 //------------------------------------------------------------------------------
 // Read exact solution
 //------------------------------------------------------------------------------
 void Configuration::read_exact (Reader &fin)
 {
     std::cout << "  Reading exact solution section\n";
 
     std::string input;
 
     fin.begin_section ("exact_soln");
 
     fin.entry("available");
     fin >> input;
     if(input == "yes")
         exact_available = true;
     else if(input == "no")
         exact_available = false;
     else
     {
         std::cout << "   Error:: unknown exact solution availability option: " << input << std::endl;
         exit (0);
     }
 
     if(exact_available)
     {
 //        fin.entry("density");
 //        fin.getline (input);
 //        exact_soln.add ("density", input);
 //
 //        fin.entry ("xvelocity");
 //        fin.getline (input);
 //        exact_soln.add ("xvelocity", input);
 //
 //        fin.entry ("yvelocity");
 //        fin.getline (input);
 //        exact_soln.add ("yvelocity", input);
 //
 //        fin.entry ("zvelocity");
 //        fin.getline (input);
 //        exact_soln.add ("zvelocity", input);
 //
 //        fin.entry ("pressure");
 //        fin.getline (input);
 //        exact_soln.add ("pressure", input);
 //
 //        fin.end_section ();
     }
     else
     {
         while(!fin.eos())
             fin >> input;
     }
 }
 
 //------------------------------------------------------------------------------
 // Read boundary conditions
 //------------------------------------------------------------------------------
 void Configuration::read_boundary_condition (Reader &fin)
 {
     std::cout << "  Reading boundary condition section\n";
     std::string input;
     fin.begin_section ("boundary_condition");
 
     fin.entry("bc");
     fin >> input;
     if(input == "mixedBC")
         bcname = MIXEDBC;
     else if(input == "periodic")
         bcname = PERIODIC;
     else if(input == "absorb")
         bcname = ABSORBINGBC;
     else if(input == "slipWall")
         bcname = SLIPWALL;
     else if(input == "doubleMachRefBC")
         bcname = DOUBLEMACHREF;
     else if(input == "FAC")
         bcname = FAC;
     else if(input == "OpenFlow")
         bcname = OPENFLOW;
     else if(input == "ScramJet")
         bcname = SCRAMJET;
     else if(input == "Nozzle")
         bcname = NOZZLE;
     else
     {
         std::cout << "   Error:: unknown boundary condition: " << input << std::endl;
         exit (0);
     }
 
     std::string boundary_names[] = {"xleft","xright","yleft","yright","zleft","zright"};
     for(int i=0; i<2*dim;++i)
     {
         fin.begin_section (boundary_names[i]);
         fin.entry("type");
         fin >> input;
         if(input=="periodic")
         {
             bcType[i] = PERIODIC;
         }
         else if(input=="dirichlet")
         {
             bcType[i] = DIRICHLET;
         }
         else if(input=="donothing")
         {
             bcType[i] = NEUMANN;
         }
         else if(input=="slipWall")
         {
             bcType[i] = SLIPWALL;
         }
         else if(input == "doubleMachRefBC")
             bcType[i] = DOUBLEMACHREF;
         else if(input == "absorb")
             bcType[i] = ABSORBINGBC;
         else if(input == "FAC")
             bcType[i] = FAC;
         else if(input == "OpenFlow")
             bcType[i] = OPENFLOW;
         else if(input == "ScramJet")
             bcType[i] = SCRAMJET;
         else if(input == "Nozzle")
             bcType[i] = NOZZLE;
         else
         {
             std::cout << "   Error: unknown boundary type " << input << std::endl;
             exit (0);
         }
         if(input == "dirichlet")
         {
             std::cout << "   Error: Dirichlet type BC not implemented yet "<< std::endl;
             exit (0);
         }
         fin.end_section();
     }
 
     fin.end_section ();
 
 
 }
 
 //------------------------------------------------------------------------------
 // Read output section
 //------------------------------------------------------------------------------
 void Configuration::read_output (Reader &fin)
 {
     std::cout << "  Reading output section\n";
 
     std::string input;
 
     fin.begin_section ("output");
 
 //    fin.entry("OOC_study");
 //    fin >> input;
 //    if(input == "yes")
 //        OOC = true;
 //    else if(input == "no")
 //        OOC = false;
 //    else
 //    {
 //        std::cout << "   Error:: unknown OOC_study option: " << input << std::endl;
 //        exit (0);
 //    }
 
 
 
 //    if(OOC && (Nx_cell.size() <2 || Ny_cell.size() < 2))
 //    {
 //        std::cout << "   Error:: Cannot perform OOC study without at least two mesh levels"<< std::endl;
 //        exit (0);
 //    }
 //
 //    if(OOC && !exact_available)
 //    {
 //        std::cout << "   Error:: Cannot perform OOC study as exact solution not avaible"<<std::endl;
 //        exit (0);
 //    }
 
 //    fin.entry ("time_stamps");
 //    fin >> t_stamps;
 //    assert (t_stamps >= 1);
 
     fin.entry ("output_path");
     fin >> OUTPUT_PATH;
 
     fin.end_section ();
 }
 
 //------------------------------------------------------------------------------
 // Read Import section
 //------------------------------------------------------------------------------
 void Configuration::read_import (Reader &fin)
 {
     std::cout << "  Reading import section\n";
 
     std::string input;
 
     fin.begin_section ("import");
 
     fin.entry ("import");
     fin >> input;
     if (input == "true") {
         import_sol = true;
     }
     else{
         import_sol = false;
     }
 
     fin.entry ("import_path");
     fin >> IMPORT_PATH;
 
     fin.end_section ();
 }
\ No newline at end of file
diff --git a/ConSol/ConSol/Base/Configuration.hpp b/ConSol/ConSol/Base/Configuration.hpp
index fa4701b..177ecac 100644
--- a/ConSol/ConSol/Base/Configuration.hpp
+++ b/ConSol/ConSol/Base/Configuration.hpp
@@ -1,187 +1,187 @@
 //
 //  Configuration.hpp
 //  ConserSolver
 //
 //  Created by Fabian Moenkeberg on 04.05.17.
 //  Copyright © 2017 FabianMoenkeberg. All rights reserved.
 //
 
 #ifndef Configuration_hpp
 #define Configuration_hpp
 
 
 //class ModelBase;
 
 #include <iostream>
 #include <vector>
 #include <ostream>
 #include <chrono>
 #include "reader.h"
 #include <map>
 
 enum evaluationMethod{PTWISE, MEANVALUE};
 enum bcName{INFLOW, OUTFLOW, SLIPWALL, EXACT, PERIODIC, DIRICHLET, NEUMANN, MIXEDBC, MIXEDBCCOMPONENT, SPECIAL, DOUBLEMACHREF, Q2I, T2I, ABSORBINGBC, Q2Q, FAC, OPENFLOW, SCRAMJET, NOZZLE};
 enum ModelName{EULER, BURGERS,SW,BUCKLEYLEVERETT, LINEARADVECTION, KPP, KPPX, KPPY};
 enum TimeIntegrationName{SSPRK2, SSPRK3, SSPRK5, FE};
 enum limiterName{NOLIMITER, EXTREMAPRES};
 enum ReconstrName{NORECONSTR, ENO, WENO};
 enum InterpBaseName{MULTIQUADS, POLYNOMIALS, GAUSSIANS, POLYHARMONICS, NOTUSED};
 enum PolyBasis{MONOMIALS, LAGRANGE};
 enum FluxName{LAXFRIEDRICHS};
 enum QuadratureName{GAUSSLOBATTO};
 enum initalCondName{SIN,EULERSOD, BUBBLE, InitKPP, InitKPPY, COSCOS, InitBURGERS, SHOCKX, SINSIN, InitISENTROPICVORTEX, InitEULERBLAST, InitVORTEXSHOCK, RP, InitDOUBLEMACHREF, InitCUBE, InitFAC, InitAIRFOIL, InitNOZZLE};
 enum sourceName{NOSOURCE};
-enum stencilChoice{STD, OTF, SYMM};
+enum stencilChoice{STD, OTF, SYMM, SYMMb};
 enum MeshType{TRI, QUAD, HYBRID, TRIQUAD, MIXED, QUADQUAD, RT, RQ};
 
 class Configuration {
 
     ModelName model;
     MeshType meshType;
     evaluationMethod evalMethod;
     bcName bcType[6];
     ReconstrName reconstr;
     InterpBaseName interpBase;
     TimeIntegrationName time_scheme;
     FluxName fluxScheme;
     limiterName limiter_name;
     bcName bcname;
     initalCondName initialCond;
     sourceName source;
     stencilChoice stencilchoice;
 
     bool reconstrOnCharact = 0;
     int order;
     int degPoly;
     double tMax;
     double cfl;
     double shapeParam;
     int maxNumWrite = 20;
     int max_iter = 50000;
     double min_residue;
     std::vector<double> tInclude;
 
     double xmin,xmax,ymin,ymax;
     int dim;
     std::vector<int> Nx_cell, Ny_cell;
 
     char* input_file;
 
     bool load_grid;
     bool import_sol;
     std::string grid_file;
     std::string mesh_type;
     std::string ic_file;
     std::string source_file;
     std::string OUTPUT_PATH;
     std::string IMPORT_PATH;
     std::string SOLN_DIR;
     bool ic_file_read;
 
     bool exact_available;
 
     int nThreads;
     std::map<std::string, double> constants;
 
 
 public:
     Configuration();
 
     Configuration(char* file);
     
     Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, int order, double tMax, double cfl);
     
     Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, int order, double tMax, std::vector<double> tInclude, double cfl);
     
     Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, int order, double tMax, double cfl);
     
     Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, int order, double tMax, std::vector<double> tInclude, double cfl);
     
     Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, bool reconstrOnCharact, int order, double tMax, double cfl);
     
     Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, bool reconstrOnCharact, int order, double tMax, std::vector<double> tInclude, double cfl);
 
     ~Configuration(){};
 
     ModelName getModel();
 
     MeshType getMeshType();
 
     TimeIntegrationName getTimeIntegration();
     
     int getOrder();
     
     evaluationMethod getEvalMethod();
     
     bcName getBc();
     
     ReconstrName getReconstr();
     
     InterpBaseName getInterpBase();
 
     limiterName getLimiter();
 
     FluxName  getFlux();
 
     sourceName getSource();
 
     initalCondName getInitialCond();
     
     std::vector<double> getTInclude();
     
     double getTmax();
     
     double getCfl();
     
     int getMaxNumWrite();
     
     bool getReconstrOnCharact();
 
     int getDim();
 
     double getXmax();
 
     double getXmin();
 
     double getYmax();
 
     double getYmin();
 
     double getShapeParam();
 
     std::vector<int>  getNx();
 
     std::vector<int>  getNy();
 
     void read ();
 
     std::string getGridFile();
 
     std::string getOutputPath();
 
     std::string getImportPath();
 
     bool getImportSol();
 
     bool getLoadGrid();
 
     int getNthreads();
 
     stencilChoice getStencilChoice();
 
     std::map<std::string, double> getConstants();
 
     int getDegPoly();
 
 private:
     void read_grid (Reader&);
     void read_constants (Reader&);
     void read_numeric (Reader&);
     void read_model (Reader&);
     void read_initial_condition (Reader&);
     void read_exact (Reader&);
     void read_boundary_condition (Reader&);
     void read_output (Reader&);
     void read_import (Reader &);
 };
 #endif /* Configuration_hpp */
diff --git a/ConSol/ConSol/Base/Setup.cpp b/ConSol/ConSol/Base/Setup.cpp
index 17d04fd..4cd4d3b 100644
--- a/ConSol/ConSol/Base/Setup.cpp
+++ b/ConSol/ConSol/Base/Setup.cpp
@@ -1,448 +1,452 @@
 //
 //  Setup.cpp
 //  ConSol
 //
 //  Created by Fabian Moenkeberg on 30.06.17.
 //  Copyright © 2017 EPFL. All rights reserved.
 //
 
 #include "Setup.hpp"
 
 Setup::Setup(Configuration config){
     dim = config.getDim();
 
     switch (config.getModel()){
         case LINEARADVECTION:
             model = new LinAdv();
             break;
         case BURGERS:
             model = new Burgers();
             break;
         case BUCKLEYLEVERETT:
             model = new BuckleyLeverett();
             break;
         case EULER:
             if (dim == 1){
                 model = new Euler();
             }else{
                 model = new Euler2D();
             }
 
             break;
         case KPP:
             model = new Kpp();
             break;
         case KPPX:
             model = new Kppx();
             break;
         case KPPY:
             model = new Kppy();
             break;
         default:
             assert(false && "Unknown Model!");
             break;
     }
 
     switch (config.getLimiter()){
         case NOLIMITER:
             limiter = new NoLimiter();
             break;
         case EXTREMAPRES:
             std::cout << "Not implemented limiter" << input_file << std::endl;
             assert(0 && "Unknown Limiter");
             break;
         default:
             assert(0 && "Unknown Limiter!");
             break;
     }
 
     switch (config.getReconstr()){
         case ENO:
             if (dim ==1){
                 reconstr = new class ENORec();
             }else if (dim == 2){
 //                reconstr = new class ENO_2D();
                 switch (config.getStencilChoice()) {
                     case STD:
                         reconstr = new class ENO_2D_std();
                         break;
                     case OTF:
                         reconstr = new class ENO_2D_otf();
                         break;
                     case SYMM:
                         reconstr = new class ENO_2D_symm();
                         break;
                 }
             }
             break;
         case WENO:
             if (dim ==1){
                 assert(0 && "Unknown Reconstruction!");
             }else if (dim == 2){
                 if (config.getMeshType() == QUAD) {
                     reconstr = new class WENO_structured();
                 }else if (config.getMeshType() == HYBRID){
                     reconstr = new class HybridRBFWENO();
-                }else{
-//                    reconstr = new class MRWENO();
+                }else if (config.getStencilChoice() == SYMM) {
+                    reconstr = new class MRWENO();
+                }else if (config.getStencilChoice() == SYMMb){
                     reconstr = new class MRWENO_b();
                 }
             }
             break;
         case NORECONSTR:
             if (dim ==1){
                 reconstr = new NoReconstr();
             }else if (dim == 2) {
                 reconstr = new NoReconstr_2D();
             }
             break;
         default:
             assert(0 && "Unknown Reconstruction!");
             break;
     }
 
     switch (config.getFlux()){
         case LAXFRIEDRICHS:
             flux = new LaxFried(config.getDim());
             break;
         default:
             assert(0 && "Unknown Flux!");
             break;
     }
 //    timeInt = new IntBase(config.getTimeIntegration(), flux);
 
     switch (config.getTimeIntegration()){
         case SSPRK2:
             timeInt = new Ssprk2(flux);
             break;
         case SSPRK3:
             timeInt = new Ssprk3(flux);
             break;
         case FE:
             timeInt = new Fe(flux);
             break;
         case SSPRK5:
             timeInt = new Ssprk5(flux);
             break;
     }
 
     switch (config.getInterpBase()){
         case MULTIQUADS:
             if (dim == 1) {
                 interpBase = new Multiquadratics(config.getOrder()-1, config.getShapeParam(), config.getDim());
             }else if (dim == 2){
 //                interpBase = new Multiquadratics_2D(config.getOrder()-1, config.getDegPoly(), config.getShapeParam());
                 if (config.getMeshType()== MIXED){
                     interpBase = new Multiquadratics_2D_mixed(config.getOrder()-1, config.getDegPoly(), config.getShapeParam());
                 }else {
                     switch (config.getStencilChoice()) {
                         case STD:
                             interpBase = new Multiquadratics_2D(config.getOrder() - 1, config.getDegPoly(),
                                                                 config.getShapeParam());
                             break;
                         case OTF:
                             interpBase = new Multiquadratics_2D_otf(config.getOrder() - 1, config.getDegPoly(),
                                                                     config.getShapeParam());
                             break;
                         case SYMM:
 //                            interpBase = new Multiquadratics_2D_otf(config.getOrder() - 1, config.getDegPoly(),
 //                                                                    config.getShapeParam());
 //                            interpBase = new Multiquadratics_2D(config.getOrder() - 1, config.getDegPoly(),
 //                                                                    config.getShapeParam());
-//                            interpBase = new Multiquadratics_2D_WENO(config.getOrder() - 1, config.getDegPoly(),
-//                                                                config.getShapeParam());
+                            interpBase = new Multiquadratics_2D_WENO(config.getOrder() - 1, config.getDegPoly(),
+                                                                config.getShapeParam());
+//                            interpBase = new Multiquadratics_2D_WENO_b(config.getOrder() - 1, config.getDegPoly(),
+//                                                                     config.getShapeParam());
+                            break;
+                        case SYMMb:
                             interpBase = new Multiquadratics_2D_WENO_b(config.getOrder() - 1, config.getDegPoly(),
                                                                      config.getShapeParam());
-                            break;
                     }
                 }
             }
             break;
         case POLYHARMONICS:
             if (dim == 1) {
                 interpBase = new Polyharmonics(config.getOrder()-1, config.getDim());
             }else if (dim == 2){
                 interpBase = new Polyharmonics_2D(config.getOrder()-1);
             }
             break;
         case GAUSSIANS:
             if (dim == 1) {
                 assert(0 && "Gaussians exist only for 2 space dimensions!");
             }else if (dim == 2){
                 interpBase = new Gaussians_2D(config.getOrder()-1);
             }
             break;
         default:
             assert(0 && "Unknown Interpolation Base!");
             break;
     }
 
     // Construct Mesh
     switch (config.getBc()){
         case PERIODIC:
             if (dim==1){
                 bc = new Periodic();}
             else
                 bc = new Periodic2();
             break;
         case ABSORBINGBC:
             if (dim==1){
                 assert(false && "Not Implemented, yet!");}
             else
                 bc = new AbsorbingBC();
             break;
         case MIXEDBC:
             if (dim==1){
                 assert(false && "Not Implemented, yet!");
                 bc = new MixedBC();}
             else
                 bc = new MixedBC_2D();
             break;
         case NEUMANN:
             assert(false && "Not Implemented, yet!");
             break;
         case SLIPWALL:
             if (dim==1){
                 assert(false && "Not Implemented, yet!");}
             else
                 bc = new SlipWallBC();
             break;
         case DOUBLEMACHREF:
             if (dim==1){
                 assert(false && "Not Implemented, yet!");}
             else
                 bc = new DoubleMachReflectionBC();
             break;
         case MIXEDBCCOMPONENT:
             assert(false && "Not Implemented, yet!");
             break;
         case DIRICHLET:
             assert(false && "Not Implemented, yet!");
             break;
         case FAC:
             if (dim==1){
                 assert(false && "Not Implemented, yet!");}
             else
                 bc = new BCFAC();
             break;
         case OPENFLOW:
             if (dim==1){
                 assert(false && "Not Implemented, yet!");}
             else
                 bc = new OpenFlowBC();
             break;
         case NOZZLE:
             if (dim==1){
                 assert(false && "Not Implemented, yet!");}
             else
                 bc = new NozzleBC();
             break;
         case SCRAMJET:
             if (dim==1){
                 assert(false && "Not Implemented, yet!");}
             else if (config.getInitialCond() == InitNOZZLE){
                 bc = new ScramJetBC(0.5);
             }
             else
                 bc = new ScramJetBC();
             break;
     }
 
     switch (config.getDim()){
         case 1:{
 
             double delta = (config.getXmax()-config.getXmin())/config.getNx()[0];
             std::vector<int> Nx = config.getNx();
             std::vector<double> x0(Nx[0]+1);
             x0.resize(Nx[0]+2);
             for (int i = 0; i < Nx[0]+2; i++){
                 x0[i] = i*delta - 0.5/Nx[0];
             }
 
             int nX0 = (int)x0.size();
             std::vector<std::vector<double>> x(nX0-1,std::vector<double>(1));
             for (int i = 0; i < nX0-1; i++){
                 x[i][0] = x0[i+1];
             }
 
             mesh = new Rect(x);
             break;}
         case 2:{
             if (!config.getLoadGrid()){
                 assert(false && "Two Dim grids must get loaded from file.");
             }
 
             switch (config.getMeshType()) {
                 case TRI:
                     mesh = new Tri();
                     break;
                 case HYBRID:
                     mesh = new FastHybridMesh();
                     break;
                 case QUAD:
                     mesh = new StructuredMesh();
                     break;
                 case MIXED:
                     mesh = new HybridMesh();
                     break;
                 default:
                     assert(0 && "Unknown Mesh Type!");
                     break;
             }
 
             break;}
     }
 
 
     switch (config.getSource()){
         case NOSOURCE:
             source = new Source();
             break;
     }
 
     switch (config.getInitialCond()){
         case SIN:
             initialCond = new Sin();
             break;
         case BUBBLE:
             initialCond = new bubble();
             break;
         case COSCOS:
             initialCond = new coscos();
             break;
         case SINSIN:
             initialCond = new sinsin();
             break;
         case InitKPP:
             initialCond = new initKPP();
             break;
         case InitKPPY:
             initialCond = new initKPPY();
             break;
         case EULERSOD:
             initialCond = new InitEulerSod();
             break;
         case InitBURGERS:
             initialCond = new initBurgers();
             break;
         case SHOCKX:
             initialCond = new shockX();
             break;
         case InitCUBE:
             initialCond = new cube();
             break;
         case InitISENTROPICVORTEX:
             initialCond = new initIsentropicVortex();
             break;
         case InitEULERBLAST:
             initialCond = new initEulerBlast();
             break;
         case InitVORTEXSHOCK:
             initialCond = new initVortexShock();
             break;
         case RP:
             initialCond = new initRP(config.getConstants());
             break;
         case InitDOUBLEMACHREF:
             initialCond = new initDoubleMachReflection();
             break;
         case InitFAC:
             if (config.getBc() == SCRAMJET){
                 initialCond = new initFAC(0, 3.0);
             }else {
                 initialCond = new initFAC();
             }
             break;
         case InitAIRFOIL:
             initialCond = new AirFoil();
             break;
         case InitNOZZLE:
             initialCond = new initNozzle(0, 0.5);
             break;
     }
 }
 
 Setup::Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Source *source, Reconstruct *reconstr, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond){
     this->mesh = mesh;
     this->model = model;
     this->timeInt = timeInt;
     this->source = source;
     this->reconstr = reconstr;
     this->flux = flux;
     this->initialCond = initialCond;
     this->limiter = limiter;
     this->interpBase = new InterpBase();
     this->bc = bc;
 }
 
 Setup::Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Source *source, Reconstruct *reconstr, InterpBase *interpBase, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond){
     this->mesh = mesh;
     this->model = model;
     this->timeInt = timeInt;
     this->source = source;
     this->reconstr = reconstr;
     this->flux = flux;
     this->initialCond = initialCond;
     this->limiter = limiter;
     this->interpBase = interpBase;
     this->bc = bc;
 }
 
 Setup::Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Reconstruct *reconstr, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond){
     this->mesh = mesh;
     this->model = model;
     this->timeInt = timeInt;
     this->source = new Source();
     this->reconstr = reconstr;
     this->flux = flux;
     this->initialCond = initialCond;
     this->limiter = limiter;
     this->interpBase = new InterpBase();
     this->bc = bc;
 }
 
 Setup::Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Reconstruct *reconstr, InterpBase *interpBase, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond){
     this->mesh = mesh;
     this->model = model;
     this->timeInt = timeInt;
     this->source = new Source();
     this->reconstr = reconstr;
     this->flux = flux;
     this->initialCond = initialCond;
     this->limiter = limiter;
     this->interpBase = interpBase;
     this->bc = bc;
 }
 
 Solution* Setup::initialize(Configuration config){
 
     model->initialize(config);
     mesh->initialize(bc, config);
     interpBase->initialize(config);
     reconstr->initialize(model, interpBase, limiter, mesh, config);
     flux->initialize(mesh, source, model, reconstr);
     limiter->initialize(interpBase, model);
     if (source->getIsUsed()){
         source->initialize(config);
     }
     timeInt->initialize(config, mesh);
     this->soln =  Solution(config, mesh, initialCond);
     
     return &soln;
 }
 
 MeshBase* Setup::getMesh(){return mesh;};
 
 BcBase* Setup::getBC() {return bc;};
 
 NumFlux* Setup::getFlux(){return flux;};
 
 ModelBase* Setup::getModel(){return model;};
 
 IntBase* Setup::getTimeInt(){return timeInt;};
 
 LimiterBase* Setup::getLimiter(){return limiter;};
 
 void Setup::exportStencils(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, Configuration config){
     reconstr->exportStencils(U, mesh, t, config.getOutputPath());
 }
diff --git a/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp
index dcac3a8..00f11e2 100644
--- a/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp
+++ b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp
@@ -1,339 +1,339 @@
 //
 // Created by Fabian Moenkeberg on 2019-12-16.
 //
 
 #include "Multiquadratics_2D_WENO.hpp"
 #include "cblas.h"
 
 Multiquadratics_2D_WENO::Multiquadratics_2D_WENO(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::Multiquadratics_2D_WENO(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::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] = xi;
     }
 
     double rad = stableRad / mesh->getApproxStencilSize(stencil);
 //    double rad = stableRad / sqrt(mesh->getSize(cellIndex));
 
     std::vector<double> ret = Multiquadratics_2D_WENO::evaluateDirect(evPoints, xB, values, 1, rad);
     return ret;
 }
 
 std::vector<double> Multiquadratics_2D_WENO::evaluateDirect(std::vector<std::vector<double>> evPoints,
                                                        std::vector<std::vector<std::vector<double>>> xB,
                                                        std::vector<double> values, double eps, double rad){
     double eps0 = 0.1/rad;
     eps = eps * rad;
     int nEv = (int) evPoints.size();
     std::vector<double> smInd(3);
     std::vector<std::vector<double>> ret(3, std::vector<double>(nEv,0));
     std::vector<double> res(nEv,0);
     int nRbf =  (int) xB.size();
     std::vector<double> coeff_old;
     int nRbf0 = std::min((int) xB.size(),3);
     int nRbf_old;
     int nP_old;
     std::vector<std::vector<std::vector<double>>> xB_;
 //    std::vector<std::vector<std::vector<double>>> xB_.reserve(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
     std::vector<double> values_;
     std::vector<double> values_old;
     int degPolyMax_ = std::max(getDegPoly(nRbf),0);
     std::vector<double> values_tmp;
     for (int degPoly_ = 0; degPoly_ <= degPolyMax_; 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;
 
         int nP2 = multi_index.size();
         nP_tmp = 0;
         for (int i = 0; i < nP; i++) {
             if (multi_index[i][0] + multi_index[i][1] <= degPoly_-2) {
                 nP_tmp++;
             }
         }
         nP2 = nP_tmp;
 //        values[2] = 14;
         if (degPoly_ == 0) {
             nRbf = 1;
             xB_.resize(nRbf,std::vector<std::vector<double>>(3, std::vector<double>(2)));
             xB_[0] = xB[0];
             values_.resize(nRbf);
             values_[0] = values[0];
         }else if (degPoly_ == 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++){
                 xB_[i] = xB[i];
                 values_[i] = values[i];
             }
 //            values_[nRbf-1] = 14;
         }else{
             nRbf = xB.size();
 //            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++){
                 xB_[i] = xB[i];
                 values_[i] = values[i];
             }
 //            values_[nRbf-1] = 14;
         }
 
         int m = nRbf + nP;
 
         int approxOrder = std::min(degPoly_ + 1, order);
 //        std::vector<double> values_tmp = values_;
         getInterpCoeff(xB_, values_, eps, nP, approxOrder);
 //        getInterpCoeff(xB_, values_tmp, eps, nP, approxOrder);
 //        getInterpCoeff(xB_, values_tmp, eps, 1, 1);
         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_] = 1.0/(std::pow(ptwiseSmInd(xB[0], xB_, values_tmp, eps),2) + eps0) + 1;
             smInd[degPoly_] = my_pow(ptwiseSmInd(xB[0], xB_, values_tmp, eps),1);
         }else{
 //            double a1,a2,a3,b1,b2,b3;
 //            if (degPolyMax_>0){
 //                for (int i = 0; i < 3; i++){
 //
 //                }
 //            }
             smInd[0] = 1;
 //            double d = std::sqrt(std::fabs(calcSize(xB[0])));
 //            double tmp = 0;
 //            double sum_tt = 0;
 //            std::vector<double>a(3);
 //            std::vector<double>tt(3);
 //            for (int jj = 1; jj <= nRbf0; jj++){
 //                a[jj-1] = (values[jj] - values[0])/d;
 //            }
 //            double tmp2 = (std::fabs(a[1]-a[0]) + std::fabs(a[2]-a[0]) + std::fabs(a[1]-a[2]))/3;
 //            for (int jj = 0; jj < nRbf0; jj++){
 //                tt[jj] = (1 + tmp2*tmp2/(1e-10 + a[jj]))/3;
 //                sum_tt += tt[jj];
 //            }
 //            for (int jj = 0; jj < nRbf0; jj++){
 //                tmp += my_pow(tt[jj]/sum_tt*a[jj],2);
 //            }
 //            tmp = tmp*std::sqrt(std::abs(calcSize(xB[0])));
 ////            std::cout << tmp << "\n";
 //            smInd[0] = tmp;
         }
         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];
         }
 
 //        if (nRbf == 6){
 //            getInterpCoeff(xB_, values, eps, nP, approxOrder);
 //            double *sVecPt = values.data();
 //            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++){
 //                for ( int j = 1; j < degPolyMax_; j++) {
 //                    ret[0][i] += ret[j][i];
 //                }
 //                ret[1][i] = y[i];
 //            }
 //        }
 
 //        values_old = values_;
         nP_old = nP;
         nRbf_old = nRbf;
 
         delete[] y;
 
         delete[] EvA;
 //        delete[] sVecPt;
     }
 //    smInd[0] = 10;
     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);
     double sum_smInd = 0;
     for ( int i = 0; i <= degPolyMax_; i++){
         sum_smInd += smInd[i];
     }
 
     for (int degPoly_ = 0; degPoly_ <= degPolyMax_; degPoly_++){
         smInd[degPolyMax_-degPoly_] /=(smInd[0]);
 //        smInd[degPolyMax_-degPoly_] = 1;
 //        if (smInd[degPolyMax_-degPoly_] < 0.1){
 //            std::cout<< smInd[degPolyMax_-degPoly_];
 //            double t = ptwiseSmInd(xB[0], xB_, values_tmp, eps);
 //        }
     }
 
     for ( int j = 0; j <= degPolyMax_; j++){
         for (int i = 0; i < nEv; i++){
 //            res[i] += (ret[degPolyMax_-j][i] - ret[degPolyMax_-j-1][i])*smInd[degPolyMax_-j]/sum_smInd*(degPolyMax_+1);
             res[i] += ret[j][i]*smInd[j];
 //            ret[degPolyMax_-j][i] += ret[degPolyMax_-j][i]*smInd[degPolyMax_-j]/sum_smInd;
         }
     }
 
 //    for (int i = 0; i < nEv; i++){
 //        res[i] += ret[0][i]*smInd[0]/sum_smInd*(degPolyMax_+1);
 //    }
 
 //    if (smInd[0] > smInd[1]&&smInd[0]> smInd[2]){
 //        res = ret[0];
 ////        std::cout << 0;
 //    }else if( smInd[1]> smInd[2]){
 //        res = ret[1];
 ////        std::cout << 1;
 //    }else{
 //        res = ret[2];
 ////        std::cout << 2;
 //    }
 
 //    std::cout << smInd[0] << ", " << smInd[1] << ", " << smInd[2] << "\n";
     return res;
 }
 
 double Multiquadratics_2D_WENO::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_+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;
             double Ds_tmp = 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 += 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(eps,j)*Ds,2) + std::pow(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::getDegPoly(int n){
     int degPoly_;
     if (false){
         degPoly_ = 1;
     }else{
         degPoly_ = (int)floor(-1.5+sqrt(1+8*(n))/2);
     }
     return degPoly_;
 }