diff --git a/HW3/CodeBlocks Fast Verifications/FFT Test/FFT Test.depend b/HW3/CodeBlocks Fast Verifications/FFT Test/FFT Test.depend index 1fa3826..3a097ee 100644 --- a/HW3/CodeBlocks Fast Verifications/FFT Test/FFT Test.depend +++ b/HW3/CodeBlocks Fast Verifications/FFT Test/FFT Test.depend @@ -1,18 +1,36 @@ # depslib dependency file v1.0 1544295319 source:/home/shaqfa/Desktop/CourseHW/SP4EHW_IM/HW3/CodeBlocks Fast Verifications/FFT Test/main.cpp "fft.h" "matrix.h" 1544295350 /home/shaqfa/Desktop/CourseHW/SP4EHW_IM/HW3/CodeBlocks Fast Verifications/FFT Test/include/fft.h "matrix.h" 1544295263 /home/shaqfa/Desktop/CourseHW/SP4EHW_IM/HW3/CodeBlocks Fast Verifications/FFT Test/include/matrix.h +1544296682 source:/home/igtomic/SP4EHW/SP4EHW_IM/HW3/CodeBlocks Fast Verifications/FFT Test/main.cpp + + "fft.h" + "matrix.h" + + +1544303585 /home/igtomic/SP4EHW/SP4EHW_IM/HW3/CodeBlocks Fast Verifications/FFT Test/include/fft.h + "matrix.h" + + + +1544296682 /home/igtomic/SP4EHW/SP4EHW_IM/HW3/CodeBlocks Fast Verifications/FFT Test/include/matrix.h + + + + + + diff --git a/HW3/CodeBlocks Fast Verifications/FFT Test/FFT Test.layout b/HW3/CodeBlocks Fast Verifications/FFT Test/FFT Test.layout index 5e8e3be..f0be894 100644 --- a/HW3/CodeBlocks Fast Verifications/FFT Test/FFT Test.layout +++ b/HW3/CodeBlocks Fast Verifications/FFT Test/FFT Test.layout @@ -1,20 +1,23 @@ - + - + - + - + + + + - + - + diff --git a/HW3/CodeBlocks Fast Verifications/FFT Test/bin/Debug/FFT Test b/HW3/CodeBlocks Fast Verifications/FFT Test/bin/Debug/FFT Test index d8fa52d..7291465 100755 Binary files a/HW3/CodeBlocks Fast Verifications/FFT Test/bin/Debug/FFT Test and b/HW3/CodeBlocks Fast Verifications/FFT Test/bin/Debug/FFT Test differ diff --git a/HW3/CodeBlocks Fast Verifications/FFT Test/include/fft.h b/HW3/CodeBlocks Fast Verifications/FFT Test/include/fft.h index baeb452..8965973 100644 --- a/HW3/CodeBlocks Fast Verifications/FFT Test/include/fft.h +++ b/HW3/CodeBlocks Fast Verifications/FFT Test/include/fft.h @@ -1,153 +1,153 @@ #ifndef FFT_HH #define FFT_HH /* ------------------------------------------------------ */ #include "matrix.h" #include /* ------------------------------------------------------ */ #include using UInt = unsigned int; using Real = double; using complex = std::complex; struct FFT { static Matrix transform(Matrix& m); static Matrix itransform(Matrix& m); static Matrix> computeFrequencies(int size); }; /* ------------------------------------------------------ */ inline Matrix FFT::transform(Matrix& m_in) { // Implementation for the Forward FFT UInt N = m_in.cols(); // The sampled size UInt Sz = N*N; // Size of the Matrix NxN fftw_plan P; // Initialize a plan // 2D Complex Matricies alocation complex *m_out_f = new complex[Sz], *m_in_f= new complex[Sz]; // Casting Values for (auto&& entry : index(m_in)){ int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); m_in_f[i*N+j] = val;} /// I've obtained this (following lines) solution from this link (Not sure but should work): /// http://www.fftw.org/fftw3_doc/Complex-numbers.html /// https://www.reddit.com/r/cpp/comments/jfk9q/is_there_a_quick_way_to_convert_fftw_complex_to_a/ P = fftw_plan_dft_2d( N, N, reinterpret_cast(m_in_f), reinterpret_cast(m_out_f), FFTW_FORWARD, FFTW_ESTIMATE // Your flags here. ); fftw_execute(P); // Casting the final output matrix -> m_out_f_f Matrix m_out_f_f(N); for (auto&& entry : index(m_out_f_f)) { int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); val = m_out_f[i*N+j];} // Dealocation (FFTW Destructor) -> Not sure about this! fftw_destroy_plan(P); fftw_free(m_in_f); fftw_free(m_out_f); // Alternative method that should de-allocate the vectors -delete [] m_in_f; -delete [] m_out_f; +//delete [] m_in_f; +//delete [] m_out_f; return m_out_f_f; } /* ------------------------------------------------------ */ inline Matrix FFT::itransform(Matrix& m_in) { // Goes the same as the previouse one // The results are not normalized in this function (Divide by N^2 the output) UInt N = m_in.cols(); UInt Sz = N*N; fftw_plan P; complex *m_out_f = new complex[Sz], *m_in_f= new complex[Sz]; for (auto&& entry : index(m_in)){ int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); m_in_f[i*N+j] = val;} P = fftw_plan_dft_2d( N, N, reinterpret_cast(m_in_f), reinterpret_cast(m_out_f), FFTW_BACKWARD, FFTW_ESTIMATE ); fftw_execute(P); Matrix m_out_f_f(N); for (auto&& entry : index(m_out_f_f)) { int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); val = m_out_f[i*N+j];} fftw_destroy_plan(P); fftw_free(m_in_f); fftw_free(m_out_f); // Alternative method that should de-allocate the vectors -delete [] m_in_f; -delete [] m_out_f; +//delete [] m_in_f; +//delete [] m_out_f; return m_out_f_f; } /* ------------------------------------------------------ */ /* ------------------------------------------------------ */ inline Matrix> FFT::computeFrequencies(int size) { /// The i stands for the real part, whereas the j stands for the imaginary part /// Explanation of procedure on the link: /// http://www.fftw.org/fftw3.pdf?fbclid=IwAR3bgzRkoSBACIJnMNfk79RX29rsVkIhWgGYg_iJcPvTM_8wLLN045ntywk // specially the figure 2.1 -- The last member is treated spatially! // x2 to run over the whole size // to account for the mirror effect of DFT symmetric matrix! (To be verified) int indicies; if (size%2 == 0) indicies = size/2; else indicies = size/2+1; // New output matrix Matrix> final_frequencies(size); for (auto&& entry : index(final_frequencies)) { int i = std::get<0>(entry); int j = std::get<1>(entry); auto& val = std::get<2>(entry); if(i < indicies) val.real(i); else val.real(i - size); if(j < indicies) val.imag(j); else val.imag(j - size); } return final_frequencies; } #endif // FFT_HH diff --git a/HW3/CodeBlocks Fast Verifications/FFT Test/obj/Debug/main.o b/HW3/CodeBlocks Fast Verifications/FFT Test/obj/Debug/main.o index b7f25d0..0686624 100644 Binary files a/HW3/CodeBlocks Fast Verifications/FFT Test/obj/Debug/main.o and b/HW3/CodeBlocks Fast Verifications/FFT Test/obj/Debug/main.o differ diff --git a/HW3/compute_temperature.cc b/HW3/compute_temperature.cc index fce0d52..1571e15 100644 --- a/HW3/compute_temperature.cc +++ b/HW3/compute_temperature.cc @@ -1,11 +1,39 @@ #include "compute_temperature.hh" #include "fft.hh" #include "material_point.hh" #include +#include "my_types.hh" /* -------------------------------------------------------------------------- */ +// Constructor for initizalizing timestep +ComputeTemperature::ComputeTemperature(Real dt) : dt(dt) {} + void ComputeTemperature::compute(System& system) { + +// For the beginning, we need to read the number of particles +UInt Sz = system.getNbParticles(); + +// Matrix is filled with the material points, depending on the system size Sz + +Matrix tempMtx(sqrt(Sz)); + for (auto&& entry : index(tempMtx)) { + int i = std::get<0>(entry); + int j = std::get<1>(entry); + auto& val = std::get<2>(entry); + Particle& par = system.getParticle(i*sqrt(Sz) + j); + auto& mp = static_cast(par); + val = mp.getTemperature(); + } + + + // FFT transform, in order to solve the equation: + // https://www.math.ubc.ca/~feldman/m267/pdeft.pdf + Matrix tempMtx_fft = FFT::transform(tempMtx); + + // Coordinates from FFT + Matrix> fft_coord = FFT::computeFrequencies(N); + } /* -------------------------------------------------------------------------- */ diff --git a/HW3/compute_temperature.hh b/HW3/compute_temperature.hh index 32f3492..0badc1b 100644 --- a/HW3/compute_temperature.hh +++ b/HW3/compute_temperature.hh @@ -1,18 +1,27 @@ #ifndef __COMPUTE_TEMPERATURE__HH__ #define __COMPUTE_TEMPERATURE__HH__ + /* -------------------------------------------------------------------------- */ #include "compute.hh" +#include "my_types.hh" -//! Compute contact interaction between ping-pong balls +//! Compute temperature class ComputeTemperature : public Compute { // Virtual implementation public: //! Penalty contact implementation void compute(System& system) override; +// Adding variables +public: + Real dt; +/*! + \dt: The stable time step for the explicit solver. +*/ + }; /* -------------------------------------------------------------------------- */ #endif //__COMPUTE_TEMPERATURE__HH__