diff --git a/source/Applications/bChannelVA/Averaging.cpp b/source/Applications/bChannelVA/Averaging.cpp
index bdc545ae7db4490baaa38f16c938905a387c2f82..20fd33178f484615f11b78189fdd6e3233c33ade 100644
--- a/source/Applications/bChannelVA/Averaging.cpp
+++ b/source/Applications/bChannelVA/Averaging.cpp
@@ -579,6 +579,7 @@ void Averaging::volumeAveragingWithMPI(double l_real)
    //#pragma omp parallel num_threads(4) //private(i)
    {
       int i = 0;
+      
 #pragma omp parallel for //private(i)//scheduler(dynamic, 1)
       for (int x3 = 0; x3 < dimensions[2]; x3++)
          for (int x2 = 0; x2 < dimensions[1]; x2++)
@@ -596,9 +597,15 @@ void Averaging::volumeAveragingWithMPI(double l_real)
                double pr = 0.0;
 
                int ll = (int)l;
+               int llz1 = ll;
+               int llz2 = ll;
+               
+               if (x3 - ll < 0) llz1 = x3;
+
+               if (x3 + ll >= dimensions[2]) llz1 = dimensions[2] - 1 - x3;
 
                //#pragma omp parallel for
-               for (int z = -ll; z <= +ll; z++)
+               for (int z = -llz1; z <= +llz2; z++)
                   for (int y = -ll; y <= +ll; y++)
                      for (int x = -ll; x <= +ll; x++)
                      {
@@ -613,7 +620,7 @@ void Averaging::volumeAveragingWithMPI(double l_real)
                         if (yy < 0)   yy = dimensions[1] + yy;
                         if (yy >= dimensions[1]) yy = yy - dimensions[1];
 
-                        if (zz < 0)   zz = 0;
+                        if (zz < 0) zz = 0;
                         if (zz >= dimensions[2]) zz = dimensions[2] - 1;
 
                         double mm = (G((double)x, l)*G((double)y, l)*G((double)z, l)) / lNorm;
@@ -942,6 +949,10 @@ void Averaging::initFluctuations()
    FlucVyMatrix.resize(dimensions[0], dimensions[1], dimensions[2], 0);
    FlucVzMatrix.resize(dimensions[0], dimensions[1], dimensions[2], 0);
    FlucPrMatrix.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaFlucVxMatrix.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaFlucVyMatrix.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaFlucVzMatrix.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaFlucPrMatrix.resize(dimensions[0], dimensions[1], dimensions[2], 0);
 }
 void Averaging::initSumOfVaFluctuations()
 {
@@ -997,7 +1008,49 @@ void Averaging::fluctuationsStress()
       XZStress[i] = vxFluc[i] * vzFluc[i];
       YZStress[i] = vyFluc[i] * vzFluc[i];
    }
-} 
+}
+void Averaging::fluctuationsStress2()
+{
+   vector<double>& vxVa = vaVxMatrix.getDataVector();
+   vector<double>& vyVa = vaVyMatrix.getDataVector();
+   vector<double>& vzVa = vaVzMatrix.getDataVector();
+   vector<double>& prVa = vaPrMatrix.getDataVector();
+
+   vector<double>& vxMean = meanVaVxMatrix.getDataVector();
+   vector<double>& vyMean = meanVaVyMatrix.getDataVector();
+   vector<double>& vzMean = meanVaVzMatrix.getDataVector();
+   vector<double>& prMean = meanVaPrMatrix.getDataVector();
+
+   vector<double>& vxFluc = vaFlucVxMatrix.getDataVector();
+   vector<double>& vyFluc = vaFlucVyMatrix.getDataVector();
+   vector<double>& vzFluc = vaFlucVzMatrix.getDataVector();
+   vector<double>& prFluc = vaFlucPrMatrix.getDataVector();
+
+   vector<double>& XXStress = vaStressXX.getDataVector();
+   vector<double>& YYStress = vaStressYY.getDataVector();
+   vector<double>& ZZStress = vaStressZZ.getDataVector();
+   vector<double>& XYStress = vaStressXY.getDataVector();
+   vector<double>& XZStress = vaStressXZ.getDataVector();
+   vector<double>& YZStress = vaStressYZ.getDataVector();
+
+   int size = (int)vxVa.size();
+
+   for (int i = 0; i < size; i++)
+   {
+      vxFluc[i] = vxVa[i] - vxMean[i];
+      vyFluc[i] = vyVa[i] - vyMean[i];
+      vzFluc[i] = vzVa[i] - vzMean[i];
+      prFluc[i] = prVa[i] - prMean[i];
+
+      XXStress[i] = vxFluc[i] * vxFluc[i];
+      YYStress[i] = vyFluc[i] * vyFluc[i];
+      ZZStress[i] = vzFluc[i] * vzFluc[i];
+      XYStress[i] = vxFluc[i] * vyFluc[i];
+      XZStress[i] = vxFluc[i] * vzFluc[i];
+      YZStress[i] = vyFluc[i] * vzFluc[i];
+   }
+}
+
 
 void Averaging::sumOfVaFluctuations()
 {
@@ -1048,10 +1101,10 @@ void Averaging::meanOfVaFluctuations(int numberOfTimeSteps)
 }
 void Averaging::writeVaFluctuationsToBinaryFiles(std::string fname, int timeStep)
 {
-   writeMatrixToBinaryFiles<double>(vaFlucVxMatrix, fname + "fluctVx" + UbSystem::toString(timeStep) + ".bin");
-   writeMatrixToBinaryFiles<double>(vaFlucVyMatrix, fname + "fluctVy" + UbSystem::toString(timeStep) + ".bin");
-   writeMatrixToBinaryFiles<double>(vaFlucVzMatrix, fname + "fluctVz" + UbSystem::toString(timeStep) + ".bin");
-   writeMatrixToBinaryFiles<double>(vaFlucPrMatrix, fname + "fluctPr" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaFlucVxMatrix, fname + "Vx" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaFlucVyMatrix, fname + "Vy" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaFlucVzMatrix, fname + "Vz" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaFlucPrMatrix, fname + "Pr" + UbSystem::toString(timeStep) + ".bin");
 }
 void Averaging::readVaFluctuationsFromBinaryFiles(std::string fname, int timeStep)
 {
@@ -1099,6 +1152,12 @@ void Averaging::initStresses()
    StressXY.resize(dimensions[0], dimensions[1], dimensions[2], 0);
    StressXZ.resize(dimensions[0], dimensions[1], dimensions[2], 0);
    StressYZ.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaStressXX.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaStressYY.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaStressZZ.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaStressXY.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaStressXZ.resize(dimensions[0], dimensions[1], dimensions[2], 0);
+   vaStressYZ.resize(dimensions[0], dimensions[1], dimensions[2], 0);
 }
 void Averaging::initSumOfVaStresses()
 {
@@ -1176,21 +1235,21 @@ void Averaging::meanOfVaStresses(int numberOfTimeSteps)
 }
 void Averaging::writeVaStressesToBinaryFiles(std::string fname, int timeStep)
 {
-   writeMatrixToBinaryFiles<double>(vaStressXX, fname + UbSystem::toString(timeStep) + "stressXX" + ".bin");
-   writeMatrixToBinaryFiles<double>(vaStressYY, fname + UbSystem::toString(timeStep) + "stressYY" + ".bin");
-   writeMatrixToBinaryFiles<double>(vaStressZZ, fname + UbSystem::toString(timeStep) + "stressZZ" + ".bin");
-   writeMatrixToBinaryFiles<double>(vaStressXY, fname + UbSystem::toString(timeStep) + "stressXY" + ".bin");
-   writeMatrixToBinaryFiles<double>(vaStressXZ, fname + UbSystem::toString(timeStep) + "stressXZ" + ".bin");
-   writeMatrixToBinaryFiles<double>(vaStressYZ, fname + UbSystem::toString(timeStep) + "stressYZ" + ".bin");
+   writeMatrixToBinaryFiles<double>(vaStressXX, fname + "XX" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaStressYY, fname + "YY" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaStressZZ, fname + "ZZ" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaStressXY, fname + "XY" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaStressXZ, fname + "XZ" + UbSystem::toString(timeStep) + ".bin");
+   writeMatrixToBinaryFiles<double>(vaStressYZ, fname + "YZ" + UbSystem::toString(timeStep) + ".bin");
 }
 void Averaging::readVaStressesFromBinaryFiles(std::string fname, int timeStep)
 {
-   readMatrixFromBinaryFiles<double>(fname + UbSystem::toString(timeStep) + "XX" + ".bin", vaStressXX);
-   readMatrixFromBinaryFiles<double>(fname + UbSystem::toString(timeStep) + "YY" + ".bin", vaStressYY);
-   readMatrixFromBinaryFiles<double>(fname + UbSystem::toString(timeStep) + "ZZ" + ".bin", vaStressZZ);
-   readMatrixFromBinaryFiles<double>(fname + UbSystem::toString(timeStep) + "XY" + ".bin", vaStressXY);
-   readMatrixFromBinaryFiles<double>(fname + UbSystem::toString(timeStep) + "XZ" + ".bin", vaStressXZ);
-   readMatrixFromBinaryFiles<double>(fname + UbSystem::toString(timeStep) + "YZ" + ".bin", vaStressYZ);
+   readMatrixFromBinaryFiles<double>(fname + "XX" + UbSystem::toString(timeStep) + ".bin", vaStressXX);
+   readMatrixFromBinaryFiles<double>(fname + "YY" + UbSystem::toString(timeStep) + ".bin", vaStressYY);
+   readMatrixFromBinaryFiles<double>(fname + "ZZ" + UbSystem::toString(timeStep) + ".bin", vaStressZZ);
+   readMatrixFromBinaryFiles<double>(fname + "XY" + UbSystem::toString(timeStep) + ".bin", vaStressXY);
+   readMatrixFromBinaryFiles<double>(fname + "XZ" + UbSystem::toString(timeStep) + ".bin", vaStressXZ);
+   readMatrixFromBinaryFiles<double>(fname + "YZ" + UbSystem::toString(timeStep) + ".bin", vaStressYZ);
 }
 void Averaging::writeMeanVaStressesToBinaryFiles(std::string fname)
 {
@@ -1211,6 +1270,12 @@ void Averaging::readMeanVaStressesFromBinaryFiles(std::string fname)
    readMatrixFromBinaryFiles<double>(fname + "YZ" + ".bin", meanVaStressYZ);
 }
 
+void Averaging::writeMeanOfVaStressesToImageFile(std::string ffname)
+{
+   array < CbArray3D<double>, 4 > matrix = { meanVaStressXX, meanVaStressYY, meanVaStressZZ, meanVaStressXY };
+   writeMatrixToImageFile(ffname, matrix);
+}
+
 //------------------------------------ planar --------------------------
 void Averaging::initPlanarAveraging()
 {
diff --git a/source/Applications/bChannelVA/Averaging.h b/source/Applications/bChannelVA/Averaging.h
index 6b963d5a433db8135593de71e8eecca931f1670a..bc968d136c611aaf002c51ce72ba5c4ed49b26ea 100644
--- a/source/Applications/bChannelVA/Averaging.h
+++ b/source/Applications/bChannelVA/Averaging.h
@@ -41,6 +41,7 @@ public:
    void initMeanOfVaFluctuations();
    void initSumOfVaFluctuations();
    void fluctuationsStress();
+   void fluctuationsStress2();
    void meanOfVaFluctuations(int numberOfTimeSteps);
    void sumOfVaFluctuations();
    void writeVaFluctuationsToBinaryFiles(std::string fname, int timeStep);
@@ -60,6 +61,7 @@ public:
    void readVaStressesFromBinaryFiles(std::string fname, int timeStep);
    void writeMeanVaStressesToBinaryFiles(std::string ffname);
    void readMeanVaStressesFromBinaryFiles(std::string ffname);
+   void writeMeanOfVaStressesToImageFile(std::string ffname);
 
    void initPlanarAveraging();
    void planarAveraging();
diff --git a/source/Applications/bChannelVA/bChannelVA.cpp b/source/Applications/bChannelVA/bChannelVA.cpp
index 7375303d9897efb1c296a1ddabcacf89a34f3979..e63b4998c6c3e1ef10e93de671ca65702a927735 100644
--- a/source/Applications/bChannelVA/bChannelVA.cpp
+++ b/source/Applications/bChannelVA/bChannelVA.cpp
@@ -10,139 +10,227 @@ using namespace std;
 //////////////////////////////////////////////////////////////////////////
 int main(int argc, char* argv[])
 {
-    try
-    {
-        SPtr<Communicator> comm = MPICommunicator::getInstance();
-        int myid = comm->getProcessID();
-
-        //Pheonix
-        //double deltaX = 1;
-        //double halfDeltaX = deltaX / 2.0;
-        //std::array<int, 3> dimensions = { 600 / (int)deltaX, 400 / (int)deltaX, 400 / (int)deltaX };
-        //std::array<double, 3> geo_origin = { halfDeltaX, halfDeltaX, halfDeltaX };
-        //std::array<double, 3> geo_spacing = { 1,1,1 };
-        //std::array<int, 6> geo_extent = { 0, dimensions[0] - 1, 0, dimensions[1] - 1, 0, dimensions[2] - 1 };
-        //double real_l = 40;
-        //double l = 40;
-
-        //int startTimeStep = 600000;
-        //int timeStep = 10000;
-        //int numberOfTimeSteps = 610000; //1200000;
-        //int numberOfSamples = numberOfTimeSteps / startTimeStep;
-        //int numberOfGridPoints = dimensions[0] * dimensions[1] * dimensions[2];
-
-
-        //Bombadil
-        string pathIn = "d:/temp/BreugemChannelAnisotrop";
-        string pathOut = "d:/temp/BreugemChannelAnisotrop";
-
-        double deltaX = 10;
-        double halfDeltaX = deltaX / 2.0;
-        std::array<int, 3> dimensions = { 600/(int)deltaX, 400/(int)deltaX, 400/(int)deltaX };
-        std::array<double, 3> geo_origin = { halfDeltaX, halfDeltaX, halfDeltaX };
-        std::array<double, 3> geo_spacing = { 10,10,10 };
-        std::array<int, 6> geo_extent = { 0, dimensions[0] - 1, 0, dimensions[1] - 1, 0, dimensions[2] - 1 };
-        double real_l = 40;
-        double l = 40;
-
-        int startTimeStep = 100;
-        int timeStep = 100;
-        int numberOfTimeSteps = 500;
-        int numberOfSamples = (numberOfTimeSteps - startTimeStep)/timeStep + 1;
-        int numberOfGridPoints = dimensions[0] * dimensions[1] * dimensions[2];
-
-        Averaging av;
-        av.setDimensions(dimensions);
-        av.setExtent(geo_extent);
-        av.setOrigin(geo_origin);
-        av.setSpacing(geo_spacing);
-        av.setDeltaX(deltaX);
-
-        av.createGeoMatrix(pathIn + "/bc/bc0.pvtu");
-        if (myid == 0) av.writeGeoMatrixToBinaryFiles(pathOut + "/va/geomatrix.bin");
-        //av.readGeoMatrixFromBinaryFiles(pathOut + "/va/geomatrix.bin");
-        av.initVolumeAveragingValues();
-
-        for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
-        {
-            av.createMQMatrix(pathIn + "/mq/mq" + UbSystem::toString(t) + ".pvtu");
-            if (myid == 0) av.writeMqMatrixToBinaryFiles(pathOut + "/va/mq/mq" + UbSystem::toString(t) + "/mq", t);
-            av.volumeAveragingWithMPI(real_l);
-            if (myid == 0)
-            {
-               av.sumOfVolumeAveragingValues();
-               av.writeVolumeAveragingValuesToBinaryFiles(pathOut + "/va/vav/vav", t);
-            }
-        }
-
-        //av.writeMqMatrixToImageFile(pathOut + "/va/vav/mq");
-        //if (myid == 0) av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", 60000);
-        //if (myid == 0) av.writeVaMatrixToImageFile(pathOut + "/va/vav/vav");
-
-        if (myid == 0)
-        {
-           av.initMeanVolumeAveragingValues();
-           av.meanOfVolumeAveragingValues(numberOfSamples);
-           av.writeMeanVolumeAveragingValuesToBinaryFiles(pathOut + "/va/mean/mean");
-           av.writeMeanMatrixToImageFile(pathOut + "/va/vav/vavMean");
-        }
-
-        //for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
-        //{
-        //   av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", t);
-        //   av.sumOfVolumeAveragingValues();
-        //}
-
-        //av.initMeanVolumeAveragingValues();
-        //av.meanOfVolumeAveragingValues(numberOfSamples);
-        //av.writeMeanVolumeAveragingValuesToBinaryFiles(pathOut + "/va/mean/mean");
-        //av.writeMeanMatrixToImageFile(pathOut + "/va/vav/vavMean");
-
-        av.initFluctuations();
-        av.initSumOfVaFluctuations();
-        av.initMeanOfVaFluctuations();
-
-        av.initStresses();
-        av.initSumOfVaStresses();
-        av.initMeanOfVaStresses();
-        av.initPlanarAveraging();
-
-        av.initVolumeAveragingFluctStressValues();
-
-        for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
-        {
-            //av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", t);
-            //av.readMqMatrixFromBinaryFiles(pathOut + "/va/mq/mq" + UbSystem::toString(t) + "/mq", t);
-            av.fluctuationsStress();
-            av.writeFluctuationsToImageFile(pathOut + "/va/vav/vaFluct");
-
-            av.volumeAveragingFluctStressWithMPI(real_l);
-            av.writeVaFluctuationsToBinaryFiles(pathOut + "/va/av/Fluc", t);
-            av.writeVaStressesToBinaryFiles(pathOut + "/va/av/Stresses", t);
-            av.writeVaFluctuationsToImageFile(pathOut + "/va/vav/vaVAFluct");
-
-            av.sumOfVaFluctuations();
-            av.sumOfVaStresses();
-        }
-
-        av.meanOfVaFluctuations(numberOfSamples);
-        av.meanOfVaStresses(numberOfSamples);
-        av.writeMeanOfVaFluctuationsToImageFile(pathOut + "/va/vav/vaMeanFluct");
-
-        av.planarAveraging();
-        av.writeToCSV(pathOut + "/va/av", geo_origin[2], deltaX);
-    }
-    catch (const std::exception& e)
-    {
-        cerr << e.what() << endl << flush;
-    }
-    catch (std::string& s)
-    {
-        cerr << s << endl;
-    }
-    catch (...)
-    {
-        cerr << "unknown exception" << endl;
-    }
+   try
+   {
+      SPtr<Communicator> comm = MPICommunicator::getInstance();
+      int myid = comm->getProcessID();
+
+      //Pheonix
+      //double deltaX = 1;
+      //double halfDeltaX = deltaX / 2.0;
+      //std::array<int, 3> dimensions = { 600 / (int)deltaX, 400 / (int)deltaX, 400 / (int)deltaX };
+      //std::array<double, 3> geo_origin = { halfDeltaX, halfDeltaX, halfDeltaX };
+      //std::array<double, 3> geo_spacing = { 1,1,1 };
+      //std::array<int, 6> geo_extent = { 0, dimensions[0] - 1, 0, dimensions[1] - 1, 0, dimensions[2] - 1 };
+      //double real_l = 40;
+      //double l = 40;
+
+      //int startTimeStep = 600000;
+      //int timeStep = 10000;
+      //int numberOfTimeSteps = 610000; //1200000;
+      //int numberOfSamples = numberOfTimeSteps / startTimeStep;
+      //int numberOfGridPoints = dimensions[0] * dimensions[1] * dimensions[2];
+
+
+      //Bombadil
+      string pathIn = "d:/temp/BreugemChannelAnisotrop2";
+      string pathOut = "d:/temp/BreugemChannelAnisotrop2";
+
+      double deltaX = 10;
+      double halfDeltaX = deltaX / 2.0;
+      std::array<int, 3> dimensions = { 600 / (int)deltaX, 400 / (int)deltaX, 400 / (int)deltaX };
+      std::array<double, 3> geo_origin = { halfDeltaX, halfDeltaX, halfDeltaX };
+      std::array<double, 3> geo_spacing = { 10,10,10 };
+      std::array<int, 6> geo_extent = { 0, dimensions[0] - 1, 0, dimensions[1] - 1, 0, dimensions[2] - 1 };
+      double real_l = 40;
+      double l = 40;
+
+      int startTimeStep = 60000;
+      int timeStep = 1000;
+      int numberOfTimeSteps = 65000;
+      int numberOfSamples = (numberOfTimeSteps - startTimeStep) / timeStep + 1;
+      int numberOfGridPoints = dimensions[0] * dimensions[1] * dimensions[2];
+
+      //Averaging av;
+      //av.setDimensions(dimensions);
+      //av.setExtent(geo_extent);
+      //av.setOrigin(geo_origin);
+      //av.setSpacing(geo_spacing);
+      //av.setDeltaX(deltaX);
+
+      //av.createGeoMatrix(pathIn + "/bc/bc0.pvtu");
+      //if (myid == 0) av.writeGeoMatrixToBinaryFiles(pathOut + "/va/geomatrix.bin");
+      ////av.readGeoMatrixFromBinaryFiles(pathOut + "/va/geomatrix.bin");
+      //av.initVolumeAveragingValues();
+
+      //for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
+      //{
+      //    av.createMQMatrix(pathIn + "/mq/mq" + UbSystem::toString(t) + ".pvtu");
+      //    if (myid == 0) av.writeMqMatrixToBinaryFiles(pathOut + "/va/mq/mq" + UbSystem::toString(t) + "/mq", t);
+      //    av.volumeAveragingWithMPI(real_l);
+      //    if (myid == 0)
+      //    {
+      //       av.sumOfVolumeAveragingValues();
+      //       av.writeVolumeAveragingValuesToBinaryFiles(pathOut + "/va/vav/vav", t);
+      //    }
+      //}
+
+      ////av.writeMqMatrixToImageFile(pathOut + "/va/vav/mq");
+      ////if (myid == 0) av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", 60000);
+      ////if (myid == 0) av.writeVaMatrixToImageFile(pathOut + "/va/vav/vav");
+
+      //if (myid == 0)
+      //{
+      //   av.initMeanVolumeAveragingValues();
+      //   av.meanOfVolumeAveragingValues(numberOfSamples);
+      //   av.writeMeanVolumeAveragingValuesToBinaryFiles(pathOut + "/va/mean/mean");
+      //   av.writeMeanMatrixToImageFile(pathOut + "/va/vav/vavMean");
+      //}
+
+      ////for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
+      ////{
+      ////   av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", t);
+      ////   av.sumOfVolumeAveragingValues();
+      ////}
+
+      ////av.initMeanVolumeAveragingValues();
+      ////av.meanOfVolumeAveragingValues(numberOfSamples);
+      ////av.writeMeanVolumeAveragingValuesToBinaryFiles(pathOut + "/va/mean/mean");
+      ////av.writeMeanMatrixToImageFile(pathOut + "/va/vav/vavMean");
+
+      //av.initFluctuations();
+      //av.initSumOfVaFluctuations();
+      //av.initMeanOfVaFluctuations();
+
+      //av.initStresses();
+      //av.initSumOfVaStresses();
+      //av.initMeanOfVaStresses();
+      //av.initPlanarAveraging();
+
+      //av.initVolumeAveragingFluctStressValues();
+
+      //for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
+      //{
+      //    //av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", t);
+      //    //av.readMqMatrixFromBinaryFiles(pathOut + "/va/mq/mq" + UbSystem::toString(t) + "/mq", t);
+      //    av.fluctuationsStress();
+      //    av.writeFluctuationsToImageFile(pathOut + "/va/vav/vaFluct");
+
+      //    av.volumeAveragingFluctStressWithMPI(real_l);
+      //    av.writeVaFluctuationsToBinaryFiles(pathOut + "/va/av/Fluc", t);
+      //    av.writeVaStressesToBinaryFiles(pathOut + "/va/av/Stresses", t);
+      //    av.writeVaFluctuationsToImageFile(pathOut + "/va/vav/vaVAFluct");
+
+      //    av.sumOfVaFluctuations();
+      //    av.sumOfVaStresses();
+      //}
+
+      //av.meanOfVaFluctuations(numberOfSamples);
+      //av.meanOfVaStresses(numberOfSamples);
+      //av.writeMeanOfVaFluctuationsToImageFile(pathOut + "/va/vav/vaMeanFluct");
+
+      //av.planarAveraging();
+      //av.writeToCSV(pathOut + "/va/av", geo_origin[2], deltaX);
+
+      Averaging av;
+
+      av.setDimensions(dimensions);
+      av.setExtent(geo_extent);
+      av.setOrigin(geo_origin);
+      av.setSpacing(geo_spacing);
+      av.setDeltaX(deltaX);
+
+      av.readGeoMatrixFromBinaryFiles(pathOut + "/va/geomatrix.bin");
+
+      av.initFluctuations();
+      av.initSumOfVaFluctuations();
+      av.initMeanOfVaFluctuations();
+
+      av.initStresses();
+      av.initSumOfVaStresses();
+      av.initMeanOfVaStresses();
+      av.initPlanarAveraging();
+
+      //av.initVolumeAveragingFluctStressValues();
+
+      av.readMeanVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/mean/mean");
+      av.writeMeanMatrixToImageFile(pathOut + "/va/vav/vavMean");
+
+      //for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
+      //{
+      //   av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", t);
+      //   av.readMqMatrixFromBinaryFiles(pathOut + "/va/mq/mq" + UbSystem::toString(t) + "/mq", t);
+      //   av.fluctuationsStress();
+      //   //av.writeFluctuationsToImageFile(pathOut + "/va/vav/vaFluct");
+
+      //   av.volumeAveragingFluctStressWithMPI(real_l);
+      //   av.writeVaFluctuationsToBinaryFiles(pathOut + "/va/vaFluc/vaFluc" + UbSystem::toString(t) + "/vaFluc", t);
+      //   av.writeVaStressesToBinaryFiles(pathOut + "/va/vaStress/vaStress" + UbSystem::toString(t) + "/vaStress", t);
+      //   //av.writeVaFluctuationsToImageFile(pathOut + "/va/vav/vaVAFluct");
+
+      //   //av.sumOfVaFluctuations();
+      //   //av.sumOfVaStresses();
+      //}
+
+      av.initVolumeAveragingValues();
+
+      for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
+      {
+         av.createMQMatrix(pathIn + "/mq/mq" + UbSystem::toString(t) + ".pvtu");
+         //av.readMqMatrixFromBinaryFiles(pathOut + "/va/mq/mq" + UbSystem::toString(t) + "/mq", t);
+         av.volumeAveragingWithMPI(real_l);
+         if (myid == 0)
+         {
+            av.sumOfVolumeAveragingValues();
+            av.writeVolumeAveragingValuesToBinaryFiles(pathOut + "/va/vav/vav", t);
+            av.writeVaMatrixToImageFile(pathOut + "/va/vav/vav" + UbSystem::toString(t));
+         }
+      }
+
+      if (myid == 0)
+      {
+         av.initMeanVolumeAveragingValues();
+         av.meanOfVolumeAveragingValues(numberOfSamples);
+         av.writeMeanVolumeAveragingValuesToBinaryFiles(pathOut + "/va/mean/mean");
+         av.writeMeanMatrixToImageFile(pathOut + "/va/vav/vavMean");
+      }
+
+      av.initVolumeAveragingFluctStressValues();
+
+      for (int t = startTimeStep; t <= numberOfTimeSteps; t += timeStep)
+      {
+         av.readVolumeAveragingValuesFromBinaryFiles(pathOut + "/va/vav/vav", t);
+         //av.readMqMatrixFromBinaryFiles(pathOut + "/va/mq/mq" + UbSystem::toString(t) + "/mq", t);
+         av.fluctuationsStress2();
+         //av.writeFluctuationsToImageFile(pathOut + "/va/vav/vaFluct");
+
+         //av.volumeAveragingFluctStressWithMPI(real_l);
+         av.writeVaFluctuationsToBinaryFiles(pathOut + "/va/vaFluc/vaFluc" + UbSystem::toString(t) + "/vaFluc", t);
+         av.writeVaStressesToBinaryFiles(pathOut + "/va/vaStress/vaStress" + UbSystem::toString(t) + "/vaStress", t);
+         //av.writeVaFluctuationsToImageFile(pathOut + "/va/vav/vaVAFluct");
+
+         av.sumOfVaFluctuations();
+         av.sumOfVaStresses();
+      }
+
+      av.meanOfVaFluctuations(numberOfSamples);
+      av.meanOfVaStresses(numberOfSamples);
+      av.writeMeanOfVaFluctuationsToImageFile(pathOut + "/va/vav/vaMeanFluct");
+
+      av.planarAveraging();
+      av.writeToCSV(pathOut + "/va/av", geo_origin[2], deltaX);
+   }
+   catch (const std::exception& e)
+   {
+      cerr << e.what() << endl << flush;
+   }
+   catch (std::string& s)
+   {
+      cerr << s << endl;
+   }
+   catch (...)
+   {
+      cerr << "unknown exception" << endl;
+   }
 
 }
\ No newline at end of file