diff --git a/hw3-heat-fft/heat_distribution.py b/hw3-heat-fft/heat_distribution.py
new file mode 100644
index 00000000..43bb6324
--- /dev/null
+++ b/hw3-heat-fft/heat_distribution.py
@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+@author: Lars Blatny and Bertil Trottet
+"""
+
+
+import numpy as np
+
+#########################################################################
+N = 11
+L = 2.0
+R = 1.0
+#########################################################################
+
+x = np.linspace(-L/2,L/2,N)
+y = np.linspace(-L/2,L/2,N)
+xv, yv = np.meshgrid(x, y, sparse=False, indexing='ij')
+hv = np.zeros((N,N))
+
+for i in range(N):
+    for j in range(N):
+        if (xv[i][j]**2+yv[i][j]**2<R**2):
+            hv[i][j] = 1.0
+
+#########################################################################
+
+f = open('hv.txt', "w")
+for i in range(N):
+    for j in range(N):
+        f.write(str(round(hv[i][j],3)))
+        f.write('\t')
+    f.write('\n')
+f.close()
+
+
+
+        
\ No newline at end of file
diff --git a/hw3-heat-fft/hv.txt b/hw3-heat-fft/hv.txt
new file mode 100644
index 00000000..64531693
--- /dev/null
+++ b/hw3-heat-fft/hv.txt
@@ -0,0 +1,11 @@
+0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	
+0.0	0.0	0.0	1.0	1.0	1.0	1.0	1.0	0.0	0.0	0.0	
+0.0	0.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.0	0.0	
+0.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.0	
+0.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.0	
+0.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.0	
+0.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.0	
+0.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.0	
+0.0	0.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	0.0	0.0	
+0.0	0.0	0.0	1.0	1.0	1.0	1.0	1.0	0.0	0.0	0.0	
+0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	
diff --git a/hw3-heat-fft/test_temperature.cc b/hw3-heat-fft/test_temperature.cc
index 45739dca..482c30ff 100644
--- a/hw3-heat-fft/test_temperature.cc
+++ b/hw3-heat-fft/test_temperature.cc
@@ -1,89 +1,130 @@
 #include "compute_temperature.hh"
 #include "csv_reader.hh"
 #include "csv_writer.hh"
 #include "material_point.hh"
 #include "material_points_factory.hh"
 #include "system.hh"
 #include <gtest/gtest.h>
 #include <vector>
 #include "vector.hh"
 
 
 
-class HomogeneousTemperatures : public ::testing::Test {};
+class HomogeneousTemperatures : public ::testing::Test {
 
-class VolumetricHeatSource : public ::testing::Test {};
+};
 
+class VolumetricHeatSource : public ::testing::Test {
+
+};
 
-TEST_F(HomogeneousTemperatures, homo_temp) {
 
+TEST_F(HomogeneousTemperatures, homo_temp) {
     System system;
     MaterialPointsFactory::getInstance();
     std::vector<MaterialPoint> materialpoints;
     UInt N = 11; // discretization of particles
     Real L = 2.0;
     Real dx = L / (N - 1.0);
     for (UInt i = 0; i < N; ++i) {
         for (UInt j = 0; j < N; ++j) {
             UInt k = i + j * N;
-            std::cout << "k = " << k << std::endl;
+            //std::cout << "k = " << k << std::endl;
             MaterialPoint p;
             Vector pos;
             pos[0] = i*dx - L/2.;
             pos[1] = j*dx - L/2.;
             pos[2] = 0;
             p.setPosition(pos);
             p.getTemperature() = 1.;
             p.getHeatRate() = 0.;
             materialpoints.push_back(p);
-            std::cout << "MaterialPoint = " << p << std::endl;
+            //std::cout << "MaterialPoint = " << p << std::endl;
             system.addParticle(std::make_shared<MaterialPoint>(p));
         } // end for j
     } // end for i
 
     auto temperature = std::make_shared<ComputeTemperature>();
     Real dt = 0.001;
     temperature->setDeltaT(dt);
     temperature->compute(system);
+
+    EXPECT_EQ(N*N, system.getNbParticles()); // CHECK THE NUMBER OF PARTICLES
+    for(auto& p : system){
+        MaterialPoint & mp = dynamic_cast<MaterialPoint&>(p);
+        Real temp = mp.getTemperature();
+        std::cout << "Temperature mp" << ":  " << temp << std::endl;
+        ASSERT_NEAR(1, temp, 1e-15); // CHECK THE TEMPERATURE VALUE FOR EACH PARTICLE
+    }
 }
 
 TEST_F(VolumetricHeatSource, volumetric_temp) {
-/*
     System system;
     MaterialPointsFactory::getInstance();
     std::vector <MaterialPoint> materialpoints;
-    n_materialpoints = 10;
-    for (UInt i = 0; i < n_materialpoints; ++i) {
-        MaterialPoint p;
-        Vector pos[2];
+    UInt N = 11; // discretization of particles
+    Real L = 2.0;
+    Real dx = L / (N - 1.0);
+    Real tolerance = 0.001;
+    for (UInt i = 0; i < N; ++i) {
+        for (UInt j = 0; j < N; j++) {
+            UInt k = i + j * N;
+            //std::cout << "k = " << k << std::endl;
+            MaterialPoint p;
+            Vector pos;
+            pos[0] = i*dx - L/2.;
+            pos[1] = j*dx - L/2.;
+            pos[2] = 0;
+            p.setPosition(pos);
+            Real x = p.getPosition()[0]; //i*dx - L/2.;
+            if (x > 0.5 - tolerance/2.0 and x < 0.5 + tolerance/2.0) {
+                p.getTemperature() = x;
+                p.getHeatRate() = 1.;
+            } else if (x > -0.5 - tolerance/2.0 and x < -0.5 + tolerance/2.0) {
+                p.getTemperature() = x;
+                p.getHeatRate() = -1.;
+            } else if (x <= -0.5 - tolerance/2.0) {
+                p.getTemperature() = -x-1.0;
+                p.getHeatRate() = 0.;
+            } else if (x >= +0.5 + tolerance/2.0){
+                p.getTemperature() = -x+1.0;
+                p.getHeatRate() = 0.;
+            } else {
+                p.getTemperature() = x;
+                p.getHeatRate() = 0.;
+            };
+            materialpoints.push_back(p);
+            //std::cout << "MaterialPoint = " << p << std::endl;
+            system.addParticle(std::make_shared<MaterialPoint>(p));
+        }
+    }
 
-        x = (n_materialpoints / i) * dx - L/2;
-        pos[0] = (n_materialpoints / i) * dx - L/2;
-        pos[1] = (i % n_materialpoints) * dx - L/2;
-        p.getPosition() = pos[2];
+    auto temperature = std::make_shared<ComputeTemperature>();
+    Real dt = 0.001;
+    temperature->setDeltaT(dt);
+    temperature->compute(system);
 
+    EXPECT_EQ(N*N, system.getNbParticles()); // CHECK THE NUMBER OF PARTICLES
+    UInt index = 0;
+    for(auto& p : system){
+        MaterialPoint & mp = dynamic_cast<MaterialPoint&>(p);
+        Real temp = mp.getTemperature();
+        std::cout << "Temperature mp" << index << " : " << temp << std::endl;
+        //ASSERT_NEAR(1, temp, 1e-15); // CHECK THE TEMPERATURE VALUE FOR EACH PARTICLE
 
+        Real x = p.getPosition()[0]; //i*dx - L/2.;
         if (x > 0.5 - tolerance/2.0 and x < 0.5 + tolerance/2.0) {
-            p.getTemperature() = x;
-            p.getHeatRate() = 1.;
+            ASSERT_NEAR(x, mp.getTemperature(), 1e-15);
         } else if (x > -0.5 - tolerance/2.0 and x < -0.5 + tolerance/2.0) {
-            p.getTemperature() = x;
-            p.getHeatRate() = -1.;
+            ASSERT_NEAR(x, mp.getTemperature(), 1e-15);
         } else if (x <= -0.5 - tolerance/2.0) {
-            p.getTemperature() = -x-1;
-            p.getHeatRate() = 0.;
+            ASSERT_NEAR(-x-1.0, mp.getTemperature(), 1e-15);
         } else if (x >= +0.5 + tolerance/2.0){
-            p.getTemperature() = -x+1;
-            p.getHeatRate() = 0.;
+            ASSERT_NEAR(-x+1.0, mp.getTemperature(), 1e-15);
         } else {
-            p.getTemperature() = x;
-            p.getHeatRate() = 0.;
-        };
-
-        materialpoints.push_back(p);
-
-        std::cout << p << std::endl;
-        system.addParticle(std::make_shared<MaterialPoint>(p));
+            ASSERT_NEAR(x, mp.getTemperature(), 1e-15);
+        }
 
-*/
+        index = index + 1;
+    }
 }