diff --git a/ConSol/ConSol/Reconstr/NoReconstr_2D.cpp b/ConSol/ConSol/Reconstr/NoReconstr_2D.cpp
index b6789b1..8d35cd8 100644
--- a/ConSol/ConSol/Reconstr/NoReconstr_2D.cpp
+++ b/ConSol/ConSol/Reconstr/NoReconstr_2D.cpp
@@ -1,53 +1,40 @@
 //
 // Created by Fabian Moenkeberg on 30.07.18.
 //
 
 #include "NoReconstr_2D.hpp"
 
 void NoReconstr_2D::reconstruct(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, std::vector<std::vector<double>> &U,
         MeshBase *mesh, double t, std::vector<std::vector<double>> &Source){
 
 
     int dimSyst = mesh->getNelem();
     int nCells = mesh->getNcells();
-//            std::vector<std::vector<std::vector<double>>> UReconstr(2,std::vector<std::vector<double>>(mesh->getNedges(),std::vector<double>(dimSyst)));
     UReconstr.resize(2,std::vector<std::vector<std::vector<double>>>(1,std::vector<std::vector<double>>(mesh->getNedges(),std::vector<double>(dimSyst))));
-
+    std::vector<double> source_loc;
 
     std::vector<int> c2e_loc(4);
     int nNodes;
-//            std::cout << "\n";
     for (int ii = 0; ii < nCells; ii++){
-//        c2e_loc.resize(mesh)
         c2e_loc = mesh->getC2E_internal(ii);
         nNodes = c2e_loc.size();
         int internal_i = mesh->getInternal(ii);
-//                std::cout << "int: " << internal_i << " for "<< ii <<",";
 
         for (int j = 0; j < nNodes; j++){
-//                    std::cout << mesh->getRight(c2e_loc[j]) << ", ";
 
             if (mesh->getRight(c2e_loc[j])== internal_i){
                 UReconstr[0][c2e_loc[j]][0] = U[internal_i];
-//                        std::cout<< "left: "<< ii << ", "<< j <<"\n";
             }else{
                 UReconstr[1][c2e_loc[j]][0] = U[internal_i];
-//                        std::cout<< "right\n";
             }
         }
-//                std::cout << "\n";
+
+        std::vector<std::vector<double>> xy = mesh->getXCoordOfCell(internal_i);
+        source_loc  = source->source(U[internal_i], t,  std::vector<double>{(xy[0][0]+xy[1][0]+xy[2][0])/3,(xy[0][1]+xy[1][1]+xy[2][1])/3});
+
+        for (int k = 0; k < dimSyst; k++) {
+            Source[internal_i][k] = source_loc[k];
+        }
     }
-//            std::cout << "URecsize: " << UReconstr[0].size() << "\n";
-//            for (int i = 0; i < UReconstr[0].size(); i++){
-//                std::cout << UReconstr[0][i][0][0] << ",";
-//            }std::cout << "\n";
-//    std::cout << "URecsize: " << UReconstr[0].size() << "\n";
-//    for (int i = 0; i < UReconstr[0].size(); i++){
-//        UReconstr[0][i][0][0] = i+1;
-//    }std::cout << "\n";
     mesh->updateEdges(UReconstr, t);
-//    std::cout << "URecsize: " << UReconstr[0].size() << "\n";
-//    for (int i = 0; i < UReconstr[0].size(); i++){
-//        std::cout << UReconstr[0][i][0][0] << ",";
-//    }std::cout << "\n";
 }
diff --git a/ConSol/ConSol/Source/AxisymmEuler_std_noSwirl.cpp b/ConSol/ConSol/Source/AxisymmEuler_std_noSwirl.cpp
index d69ffb9..45ebfa5 100644
--- a/ConSol/ConSol/Source/AxisymmEuler_std_noSwirl.cpp
+++ b/ConSol/ConSol/Source/AxisymmEuler_std_noSwirl.cpp
@@ -1,253 +1,274 @@
 //
 // Created by Fabian Moenkeberg on 2020-02-13.
 //
 
 #include "AxisymmEuler_std_noSwirl.hpp"
 
 std::vector<double> AxisymmEuler_std_noSwirl::source(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, int iCell){
     std::vector<double> ret(U[iCell].size(),0);
 
     long double rho, p, m1, m2, E;
     rho = U[iCell][0];
     m1 = U[iCell][1];
     m2 = U[iCell][2];
     E = U[iCell][3];
     p = (gamma - 1)*(E - 0.5*(m1*m1 + m2*m2)/rho);
 
 //    if (iCell == 2225){
 //        std::cout << iCell;
 //    }
     std::vector<double> r = mesh->MidPoint(iCell);
     ret[0]  = -m2/r[1];
     ret[1]  = -m1*m2/rho/r[1];
     ret[2]  = -m2*m2/rho/r[1];
     ret[3]  = -m2/rho*(E+p)/r[1];
     return ret;
 }
 
 void AxisymmEuler_std_noSwirl::source(std::vector<std::vector<double>> &U, MeshBase *mesh, double t,
        std::vector<std::vector<double>> &numFlux, std::vector<std::vector<double>> &HighOrdersource){
 
     omp_set_num_threads(nThreads);
 #pragma omp parallel for
     for (int i = 0; i < nCells;i++) {
         int iCell = mesh->getInternal(i);
 //        long double rho, p, m1, m2, E;
 //        rho = U[iCell][0];
 //        m1 = U[iCell][1];
 //        m2 = U[iCell][2];
 //        E = U[iCell][3];
 //        p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
 
 //    if (iCell == 2225){
 //        std::cout << iCell;
 //    }
 //        std::vector<std::vector<double>>  xCoord = mesh->getXCoordOfCell(iCell);
 //        double r = (xCoord[0][1] + xCoord[1][1] + xCoord[2][1])/3;
 //        numFlux[iCell][0] -= m2 / r;
 //        numFlux[iCell][1] -= m1 * m2 / rho / r;
 //        numFlux[iCell][2] -= m2 * m2 / rho / r;
 //        numFlux[iCell][3] -= m2 / rho * (E + p) / r;
 
         numFlux[iCell][0] -= HighOrdersource[iCell][0];
         numFlux[iCell][1] -= HighOrdersource[iCell][1];
         numFlux[iCell][2] -= HighOrdersource[iCell][2];
         numFlux[iCell][3] -= HighOrdersource[iCell][3];
     }
 
     omp_set_num_threads(nThreads);
-//#pragma omp parallel for
-//    for (int ii = 0; ii < nDomainsQUAD; ii++){
-//        for (int i = 0; i < nX[ii]; i++) {
-//            for (int j = 0; j < nY[ii]; j++) {
-//                long double rho, p, m1, m2, E;
-//                int iCell = i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT];
-//
-//                rho = U[iCell][0];
-//                m1 = U[iCell][1];
-//                m2 = U[iCell][2];
-//                E = U[iCell][3];
-//                p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
-//
-//                double r = yMin[ii] + (j + ngc + 0.5)*dY[ii];
-//                numFlux[iCell][0] -= m2 / r;
-//                numFlux[iCell][1] -= m1 * m2 / rho / r;
-//                numFlux[iCell][2] -= m2 * m2 / rho / r;
-//                numFlux[iCell][3] -= m2 / rho * (E + p) / r;
-//            }
-//        }
-//    }
+    if (order_scheme == 1) {
 #pragma omp parallel for
-    for (int ii = 0; ii < nDomainsQUAD; ii++){
-        for (int i = 0; i < nX[ii]; i++) {
-            for (int j = 0; j < nY[ii]; j++) {
-                long double rho, p, m1, m2, E, r, w;
-                int i0 =  + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT];
-                int iCell = i+ngc+(j+ngc)*(nX[ii]+2*ngc);
-                int i00 = i+ngc-1+(j+ngc-1)*(nX[ii]+2*ngc);
-                int i01 = i+ngc-1+(j+ngc)*(nX[ii]+2*ngc);
-                int i02 = i+ngc-1+(j+ngc+1)*(nX[ii]+2*ngc);
-                int i10 = i+ngc+(j+ngc-1)*(nX[ii]+2*ngc);
-                int i11 = i+ngc+(j+ngc)*(nX[ii]+2*ngc);
-                int i12 = i+ngc+(j+ngc+1)*(nX[ii]+2*ngc);
-                int i20 = i+ngc+1+(j+ngc-1)*(nX[ii]+2*ngc);
-                int i21 = i+ngc+1+(j+ngc)*(nX[ii]+2*ngc);
-                int i22 = i+ngc+1+(j+ngc+1)*(nX[ii]+2*ngc);
-
-                // For calculating u00
-                int is = 0;
-                rho = (U[i00][is] + 4*sqrt(3)*U[i01][is] - U[i02][is] + 4*sqrt(3)*U[i10][is] + 48*U[i11][is] - 4*sqrt(3)*U[i12][is]
-                        - U[i20][is] - 4*sqrt(3)*U[i21][is] + U[i22][is])/48;
-                is = 1;
-                m1 = (U[i00][is] + 4*sqrt(3)*U[i01][is] - U[i02][is] + 4*sqrt(3)*U[i10][is] + 48*U[i11][is] - 4*sqrt(3)*U[i12][is]
-                        - U[i20][is] - 4*sqrt(3)*U[i21][is] + U[i22][is])/48;
-                is = 2;
-                m2 = (U[i00][is] + 4*sqrt(3)*U[i01][is] - U[i02][is] + 4*sqrt(3)*U[i10][is] + 48*U[i11][is] - 4*sqrt(3)*U[i12][is]
-                        - U[i20][is] - 4*sqrt(3)*U[i21][is] + U[i22][is])/48;
-                is = 3;
-                E = (U[i00][is] + 4*sqrt(3)*U[i01][is] - U[i02][is] + 4*sqrt(3)*U[i10][is] + 48*U[i11][is] - 4*sqrt(3)*U[i12][is]
-                       - U[i20][is] - 4*sqrt(3)*U[i21][is] + U[i22][is])/48;
-                p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
+        for (int ii = 0; ii < nDomainsQUAD; ii++) {
+            for (int i = 0; i < nX[ii]; i++) {
+                for (int j = 0; j < nY[ii]; j++) {
+                    long double rho, p, m1, m2, E;
+                    int iCell = i + ngc + (j + ngc) * (nX[ii] + 2 * ngc) +
+                                submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT];
 
-                r = yMin[ii] + (j + ngc + 0.5 - 1.0/sqrt(12))*dY[ii];
-                w = 0.25;
-                numFlux[iCell][0] -= w*m2 / r;
-                numFlux[iCell][1] -= w*m1 * m2 / rho / r;
-                numFlux[iCell][2] -= w*m2 * m2 / rho / r;
-                numFlux[iCell][3] -= w*m2 / rho * (E + p) / r;
-
-                // For calculating u02
-                is = 0;
-                rho = (-U[i00][is] + 4*sqrt(3)*U[i01][is] + U[i02][is] - 4*sqrt(3)*U[i10][is] + 48*U[i11][is] + 4*sqrt(3)*U[i12][is]
-                         + U[i20][is] - 4*sqrt(3)*U[i21][is] - U[i22][is])/48;
-                is = 1;
-                m1 = (-U[i00][is] + 4*sqrt(3)*U[i01][is] + U[i02][is] - 4*sqrt(3)*U[i10][is] + 48*U[i11][is] + 4*sqrt(3)*U[i12][is]
-                        + U[i20][is] - 4*sqrt(3)*U[i21][is] - U[i22][is])/48;
-                is = 2;
-                m2 = (-U[i00][is] + 4*sqrt(3)*U[i01][is] + U[i02][is] - 4*sqrt(3)*U[i10][is] + 48*U[i11][is] + 4*sqrt(3)*U[i12][is]
-                        + U[i20][is] - 4*sqrt(3)*U[i21][is] - U[i22][is])/48;
-                is = 3;
-                E = (-U[i00][is] + 4*sqrt(3)*U[i01][is] + U[i02][is] - 4*sqrt(3)*U[i10][is] + 48*U[i11][is] + 4*sqrt(3)*U[i12][is]
-                       + U[i20][is] - 4*sqrt(3)*U[i21][is] - U[i22][is])/48;
-                p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
+                    rho = U[iCell][0];
+                    m1 = U[iCell][1];
+                    m2 = U[iCell][2];
+                    E = U[iCell][3];
+                    p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
 
-                r = yMin[ii] + (j + ngc + 0.5 + 1.0/sqrt(12))*dY[ii];
-                numFlux[iCell][0] -= w*m2 / r;
-                numFlux[iCell][1] -= w*m1 * m2 / rho / r;
-                numFlux[iCell][2] -= w*m2 * m2 / rho / r;
-                numFlux[iCell][3] -= w*m2 / rho * (E + p) / r;
-
-                // For calculating u20
-                is = 0;
-                rho = (-U[i00][is] - 4*sqrt(3)*U[i01][is] + U[i02][is] + 4*sqrt(3)*U[i10][is] + 48*U[i11][is] - 4*sqrt(3)*U[i12][is]
-                       + U[i20][is] + 4*sqrt(3)*U[i21][is] - U[i22][is])/48;
-                is = 1;
-                m1 = (-U[i00][is] - 4*sqrt(3)*U[i01][is] + U[i02][is] + 4*sqrt(3)*U[i10][is] + 48*U[i11][is] - 4*sqrt(3)*U[i12][is]
-                      + U[i20][is] + 4*sqrt(3)*U[i21][is] - U[i22][is])/48;
-                is = 2;
-                m2 = (-U[i00][is] - 4*sqrt(3)*U[i01][is] + U[i02][is] + 4*sqrt(3)*U[i10][is] + 48*U[i11][is] - 4*sqrt(3)*U[i12][is]
-                      + U[i20][is] + 4*sqrt(3)*U[i21][is] - U[i22][is])/48;
-                is = 3;
-                E = (-U[i00][is] - 4*sqrt(3)*U[i01][is] + U[i02][is] + 4*sqrt(3)*U[i10][is] + 48*U[i11][is] - 4*sqrt(3)*U[i12][is]
-                     + U[i20][is] + 4*sqrt(3)*U[i21][is] - U[i22][is])/48;
-                p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
+                    double r = yMin[ii] + (j + ngc + 0.5) * dY[ii];
+                    numFlux[iCell][0] -= m2 / r;
+                    numFlux[iCell][1] -= m1 * m2 / rho / r;
+                    numFlux[iCell][2] -= m2 * m2 / rho / r;
+                    numFlux[iCell][3] -= m2 / rho * (E + p) / r;
+                }
+            }
+        }
+    }
+    else {
+#pragma omp parallel for
+        for (int ii = 0; ii < nDomainsQUAD; ii++) {
+            for (int i = 0; i < nX[ii]; i++) {
+                for (int j = 0; j < nY[ii]; j++) {
+                    long double rho, p, m1, m2, E, r, w;
+                    int i0 = +submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT];
+                    int iCell = i + ngc + (j + ngc) * (nX[ii] + 2 * ngc);
+                    int i00 = i + ngc - 1 + (j + ngc - 1) * (nX[ii] + 2 * ngc);
+                    int i01 = i + ngc - 1 + (j + ngc) * (nX[ii] + 2 * ngc);
+                    int i02 = i + ngc - 1 + (j + ngc + 1) * (nX[ii] + 2 * ngc);
+                    int i10 = i + ngc + (j + ngc - 1) * (nX[ii] + 2 * ngc);
+                    int i11 = i + ngc + (j + ngc) * (nX[ii] + 2 * ngc);
+                    int i12 = i + ngc + (j + ngc + 1) * (nX[ii] + 2 * ngc);
+                    int i20 = i + ngc + 1 + (j + ngc - 1) * (nX[ii] + 2 * ngc);
+                    int i21 = i + ngc + 1 + (j + ngc) * (nX[ii] + 2 * ngc);
+                    int i22 = i + ngc + 1 + (j + ngc + 1) * (nX[ii] + 2 * ngc);
 
-                r = yMin[ii] + (j + ngc + 0.5 - 1.0/sqrt(12))*dY[ii];
-                numFlux[iCell][0] -= w*m2 / r;
-                numFlux[iCell][1] -= w*m1 * m2 / rho / r;
-                numFlux[iCell][2] -= w*m2 * m2 / rho / r;
-                numFlux[iCell][3] -= w*m2 / rho * (E + p) / r;
-
-                // For calculating u22
-                is = 0;
-                rho = (U[i00][is] - 4*sqrt(3)*U[i01][is] - U[i02][is] - 4*sqrt(3)*U[i10][is] + 48*U[i11][is] + 4*sqrt(3)*U[i12][is]
-                       - U[i20][is] + 4*sqrt(3)*U[i21][is] + U[i22][is])/48;
-                is = 1;
-                m1 = (U[i00][is] - 4*sqrt(3)*U[i01][is] - U[i02][is] - 4*sqrt(3)*U[i10][is] + 48*U[i11][is] + 4*sqrt(3)*U[i12][is]
-                      - U[i20][is] + 4*sqrt(3)*U[i21][is] + U[i22][is])/48;
-                is = 2;
-                m2 = (U[i00][is] - 4*sqrt(3)*U[i01][is] - U[i02][is] - 4*sqrt(3)*U[i10][is] + 48*U[i11][is] + 4*sqrt(3)*U[i12][is]
-                      - U[i20][is] + 4*sqrt(3)*U[i21][is] + U[i22][is])/48;
-                is = 3;
-                E = (U[i00][is] - 4*sqrt(3)*U[i01][is] - U[i02][is] - 4*sqrt(3)*U[i10][is] + 48*U[i11][is] + 4*sqrt(3)*U[i12][is]
-                     - U[i20][is] + 4*sqrt(3)*U[i21][is] + U[i22][is])/48;
-                p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
+                    // For calculating u00
+                    int is = 0;
+                    rho = (U[i00][is] + 4 * sqrt(3) * U[i01][is] - U[i02][is] + 4 * sqrt(3) * U[i10][is] +
+                           48 * U[i11][is] - 4 * sqrt(3) * U[i12][is]
+                           - U[i20][is] - 4 * sqrt(3) * U[i21][is] + U[i22][is]) / 48;
+                    is = 1;
+                    m1 = (U[i00][is] + 4 * sqrt(3) * U[i01][is] - U[i02][is] + 4 * sqrt(3) * U[i10][is] +
+                          48 * U[i11][is] - 4 * sqrt(3) * U[i12][is]
+                          - U[i20][is] - 4 * sqrt(3) * U[i21][is] + U[i22][is]) / 48;
+                    is = 2;
+                    m2 = (U[i00][is] + 4 * sqrt(3) * U[i01][is] - U[i02][is] + 4 * sqrt(3) * U[i10][is] +
+                          48 * U[i11][is] - 4 * sqrt(3) * U[i12][is]
+                          - U[i20][is] - 4 * sqrt(3) * U[i21][is] + U[i22][is]) / 48;
+                    is = 3;
+                    E = (U[i00][is] + 4 * sqrt(3) * U[i01][is] - U[i02][is] + 4 * sqrt(3) * U[i10][is] +
+                         48 * U[i11][is] - 4 * sqrt(3) * U[i12][is]
+                         - U[i20][is] - 4 * sqrt(3) * U[i21][is] + U[i22][is]) / 48;
+                    p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
+
+                    r = yMin[ii] + (j + ngc + 0.5 - 1.0 / sqrt(12)) * dY[ii];
+                    w = 0.25;
+                    numFlux[iCell][0] -= w * m2 / r;
+                    numFlux[iCell][1] -= w * m1 * m2 / rho / r;
+                    numFlux[iCell][2] -= w * m2 * m2 / rho / r;
+                    numFlux[iCell][3] -= w * m2 / rho * (E + p) / r;
+
+                    // For calculating u02
+                    is = 0;
+                    rho = (-U[i00][is] + 4 * sqrt(3) * U[i01][is] + U[i02][is] - 4 * sqrt(3) * U[i10][is] +
+                           48 * U[i11][is] + 4 * sqrt(3) * U[i12][is]
+                           + U[i20][is] - 4 * sqrt(3) * U[i21][is] - U[i22][is]) / 48;
+                    is = 1;
+                    m1 = (-U[i00][is] + 4 * sqrt(3) * U[i01][is] + U[i02][is] - 4 * sqrt(3) * U[i10][is] +
+                          48 * U[i11][is] + 4 * sqrt(3) * U[i12][is]
+                          + U[i20][is] - 4 * sqrt(3) * U[i21][is] - U[i22][is]) / 48;
+                    is = 2;
+                    m2 = (-U[i00][is] + 4 * sqrt(3) * U[i01][is] + U[i02][is] - 4 * sqrt(3) * U[i10][is] +
+                          48 * U[i11][is] + 4 * sqrt(3) * U[i12][is]
+                          + U[i20][is] - 4 * sqrt(3) * U[i21][is] - U[i22][is]) / 48;
+                    is = 3;
+                    E = (-U[i00][is] + 4 * sqrt(3) * U[i01][is] + U[i02][is] - 4 * sqrt(3) * U[i10][is] +
+                         48 * U[i11][is] + 4 * sqrt(3) * U[i12][is]
+                         + U[i20][is] - 4 * sqrt(3) * U[i21][is] - U[i22][is]) / 48;
+                    p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
+
+                    r = yMin[ii] + (j + ngc + 0.5 + 1.0 / sqrt(12)) * dY[ii];
+                    numFlux[iCell][0] -= w * m2 / r;
+                    numFlux[iCell][1] -= w * m1 * m2 / rho / r;
+                    numFlux[iCell][2] -= w * m2 * m2 / rho / r;
+                    numFlux[iCell][3] -= w * m2 / rho * (E + p) / r;
+
+                    // For calculating u20
+                    is = 0;
+                    rho = (-U[i00][is] - 4 * sqrt(3) * U[i01][is] + U[i02][is] + 4 * sqrt(3) * U[i10][is] +
+                           48 * U[i11][is] - 4 * sqrt(3) * U[i12][is]
+                           + U[i20][is] + 4 * sqrt(3) * U[i21][is] - U[i22][is]) / 48;
+                    is = 1;
+                    m1 = (-U[i00][is] - 4 * sqrt(3) * U[i01][is] + U[i02][is] + 4 * sqrt(3) * U[i10][is] +
+                          48 * U[i11][is] - 4 * sqrt(3) * U[i12][is]
+                          + U[i20][is] + 4 * sqrt(3) * U[i21][is] - U[i22][is]) / 48;
+                    is = 2;
+                    m2 = (-U[i00][is] - 4 * sqrt(3) * U[i01][is] + U[i02][is] + 4 * sqrt(3) * U[i10][is] +
+                          48 * U[i11][is] - 4 * sqrt(3) * U[i12][is]
+                          + U[i20][is] + 4 * sqrt(3) * U[i21][is] - U[i22][is]) / 48;
+                    is = 3;
+                    E = (-U[i00][is] - 4 * sqrt(3) * U[i01][is] + U[i02][is] + 4 * sqrt(3) * U[i10][is] +
+                         48 * U[i11][is] - 4 * sqrt(3) * U[i12][is]
+                         + U[i20][is] + 4 * sqrt(3) * U[i21][is] - U[i22][is]) / 48;
+                    p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
+
+                    r = yMin[ii] + (j + ngc + 0.5 - 1.0 / sqrt(12)) * dY[ii];
+                    numFlux[iCell][0] -= w * m2 / r;
+                    numFlux[iCell][1] -= w * m1 * m2 / rho / r;
+                    numFlux[iCell][2] -= w * m2 * m2 / rho / r;
+                    numFlux[iCell][3] -= w * m2 / rho * (E + p) / r;
+
+                    // For calculating u22
+                    is = 0;
+                    rho = (U[i00][is] - 4 * sqrt(3) * U[i01][is] - U[i02][is] - 4 * sqrt(3) * U[i10][is] +
+                           48 * U[i11][is] + 4 * sqrt(3) * U[i12][is]
+                           - U[i20][is] + 4 * sqrt(3) * U[i21][is] + U[i22][is]) / 48;
+                    is = 1;
+                    m1 = (U[i00][is] - 4 * sqrt(3) * U[i01][is] - U[i02][is] - 4 * sqrt(3) * U[i10][is] +
+                          48 * U[i11][is] + 4 * sqrt(3) * U[i12][is]
+                          - U[i20][is] + 4 * sqrt(3) * U[i21][is] + U[i22][is]) / 48;
+                    is = 2;
+                    m2 = (U[i00][is] - 4 * sqrt(3) * U[i01][is] - U[i02][is] - 4 * sqrt(3) * U[i10][is] +
+                          48 * U[i11][is] + 4 * sqrt(3) * U[i12][is]
+                          - U[i20][is] + 4 * sqrt(3) * U[i21][is] + U[i22][is]) / 48;
+                    is = 3;
+                    E = (U[i00][is] - 4 * sqrt(3) * U[i01][is] - U[i02][is] - 4 * sqrt(3) * U[i10][is] +
+                         48 * U[i11][is] + 4 * sqrt(3) * U[i12][is]
+                         - U[i20][is] + 4 * sqrt(3) * U[i21][is] + U[i22][is]) / 48;
+                    p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
 
-                r = yMin[ii] + (j + ngc + 0.5 + 1.0/sqrt(12))*dY[ii];
-                numFlux[iCell][0] -= w*m2 / r;
-                numFlux[iCell][1] -= w*m1 * m2 / rho / r;
-                numFlux[iCell][2] -= w*m2 * m2 / rho / r;
-                numFlux[iCell][3] -= w*m2 / rho * (E + p) / r;
+                    r = yMin[ii] + (j + ngc + 0.5 + 1.0 / sqrt(12)) * dY[ii];
+                    numFlux[iCell][0] -= w * m2 / r;
+                    numFlux[iCell][1] -= w * m1 * m2 / rho / r;
+                    numFlux[iCell][2] -= w * m2 * m2 / rho / r;
+                    numFlux[iCell][3] -= w * m2 / rho * (E + p) / r;
+                }
             }
         }
     }
 }
 
 void AxisymmEuler_std_noSwirl::source(std::vector<std::vector<double>> &U, MeshBase *mesh, double t,
                                       std::vector<std::vector<double>> &numFlux){
 
     omp_set_num_threads(nThreads);
 #pragma omp parallel for
     for (int i = 0; i < nCells;i++) {
         long double rho, p, m1, m2, E;
         int iCell = mesh->getInternal(i);
         rho = U[iCell][0];
         m1 = U[iCell][1];
         m2 = U[iCell][2];
         E = U[iCell][3];
         p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
 
 //    if (iCell == 2225){
 //        std::cout << iCell;
 //    }
         std::vector<std::vector<double>>  xCoord = mesh->getXCoordOfCell(iCell);
         double r = (xCoord[0][1] + xCoord[1][1] + xCoord[2][1])/3;
         numFlux[iCell][0] -= m2 / r;
         numFlux[iCell][1] -= m1 * m2 / rho / r;
         numFlux[iCell][2] -= m2 * m2 / rho / r;
         numFlux[iCell][3] -= m2 / rho * (E + p) / r;
 
     }
 
     omp_set_num_threads(nThreads);
 #pragma omp parallel for
     for (int ii = 0; ii < nDomainsQUAD; ii++){
         for (int i = 0; i < nX[ii]; i++) {
             for (int j = 0; j < nY[ii]; j++) {
                 long double rho, p, m1, m2, E;
                 int iCell = i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT];
 
                 rho = U[iCell][0];
                 m1 = U[iCell][1];
                 m2 = U[iCell][2];
                 E = U[iCell][3];
                 p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
 
                 double r = yMin[ii] + (j + ngc + 0.5)*dY[ii];
                 numFlux[iCell][0] -= m2 / r;
                 numFlux[iCell][1] -= m1 * m2 / rho / r;
                 numFlux[iCell][2] -= m2 * m2 / rho / r;
                 numFlux[iCell][3] -= m2 / rho * (E + p) / r;
             }
         }
     }
 
 }
 
 std::vector<double> AxisymmEuler_std_noSwirl::source(std::vector<double> &U, double t, std::vector<double> xCoord){
     std::vector<double> ret(4);
     long double rho, p, m1, m2, E;
     rho = U[0];
     m1 = U[1];
     m2 = U[2];
     E = U[3];
     p = (gamma - 1) * (E - 0.5 * (m1 * m1 + m2 * m2) / rho);
 
     double r = xCoord[1];
     ret[0] = m2 / r;
     ret[1] = m1 * m2 / rho / r;
     ret[2] = m2 * m2 / rho / r;
     ret[3] = m2 / rho * (E + p) / r;
 
     return ret;
 }
\ No newline at end of file
diff --git a/ConSol/ConSol/Source/Source.cpp b/ConSol/ConSol/Source/Source.cpp
index 7dfae32..c7fb094 100644
--- a/ConSol/ConSol/Source/Source.cpp
+++ b/ConSol/ConSol/Source/Source.cpp
@@ -1,56 +1,57 @@
 //
 //  SourceBase.cpp
 //  ConserSolver
 //
 //  Created by Fabian Moenkeberg on 10.05.17.
 //  Copyright © 2017 FabianMoenkeberg. All rights reserved.
 //
 
 #include "Source.hpp"
 #include "../Mesh/MeshBase.hpp"
 
 Source::Source(){};
 
 Source::Source(ModelBase model){
     this->model = model;
 };
 
 void Source::initialize(Configuration config, MeshBase* mesh){
     nCells = mesh->getNcells();
     nDomainsQUAD = mesh->getNdomainsQUAD() + mesh->getNdomainsQUADQUAD() + mesh->getNdomainsRQ();
 //    nDomainsQUADQUAD = mesh->getNdomainsQUADQUAD();
 //    nDomainsRQ = mesh->getNdomainsRQ();
     nDomainsTRIQUAD = mesh->getNdomainsTRIQUAD();
     nDomainsTRI = mesh->getNdomainsTRI();
     nDomainsRT = mesh->getNdomainsRT();
     submeshCell0 = mesh->getSubmeshCell0();
     submeshEdg0 = mesh->getSubmeshEdg0();
     nX = mesh->getNx();
     nY = mesh->getNy();
     ngc = mesh->getNgc();
     xMin = mesh->getXmin();
     yMin = mesh->getYmin();
     dX = mesh->getDx();
     dY = mesh->getDy();
+    order_scheme = config.getOrder();
 
     this->nThreads = config.getNthreads();
     if (nThreads == 0){
         nThreads = omp_get_max_threads();
     }
 };
 
 std::vector<double> Source::source(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, int iCell){
     std::vector<double> ret(0);
     return ret;
 }
 
 void Source::source(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, std::vector<std::vector<double>> &numFlux,
         std::vector<std::vector<double>> &HighOrdersource){
 }
 
 std::vector<double> Source::source(std::vector<double> &U, double t, std::vector<double> xCoord){
     std::vector<double> ret(0);
     return ret;
 }
 
 bool Source::getIsUsed(){return isUsed;};
diff --git a/ConSol/ConSol/Source/Source.hpp b/ConSol/ConSol/Source/Source.hpp
index bd953b5..f2f3b28 100644
--- a/ConSol/ConSol/Source/Source.hpp
+++ b/ConSol/ConSol/Source/Source.hpp
@@ -1,59 +1,61 @@
 //
 //  SourceBase.hpp
 //  ConserSolver
 //
 //  Created by Fabian Moenkeberg on 10.05.17.
 //  Copyright © 2017 FabianMoenkeberg. All rights reserved.
 //
 
 #ifndef Source_hpp
 #define Source_hpp
 
 class MeshBase;
 
 #include <stdio.h>
 #include <vector>
 #include "../Model/ModelBase.hpp"
 #include "../Base/Configuration.hpp"
 #include <omp.h>
 
 class Source{
 protected:
     ModelBase model;
     bool isUsed = 0;
     std::vector<int> submeshCell0;
     std::vector<int> submeshEdg0;
     int nCells;
     int nDomainsQUAD;
 //    int nDomainsQUADQUAD;
 //    int nDomainsRQ;
     int nDomainsTRI;
     int nDomainsTRIQUAD;
     int nDomainsRT;
     std::vector<int> nY;
     std::vector<int> nX;
     int ngc;
     std::vector<double> dX;
     std::vector<double> dY;
     std::vector<double> xMin;
     std::vector<double> yMin;
     int nThreads;
+    int order_scheme;
+    
 public:
     Source();
 
     virtual ~Source(){};
     
     Source(ModelBase model);
     
     virtual void initialize(Configuration config, MeshBase* mesh);
     
     virtual std::vector<double> source(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, int iCell);
 
     virtual void source(std::vector<std::vector<double>> &U, MeshBase *mesh, double t, std::vector<std::vector<double>> &numFlux,
             std::vector<std::vector<double>> &HighOrdersource);
 
     virtual std::vector<double> source(std::vector<double> &U, double t, std::vector<double> xCoord);
 
     bool getIsUsed();
 };
 #endif /* Source_hpp */