diff --git a/Homework3/src/compute_temperature.cc b/Homework3/src/compute_temperature.cc index 37627a0..91ccae5 100644 --- a/Homework3/src/compute_temperature.cc +++ b/Homework3/src/compute_temperature.cc @@ -1,65 +1,65 @@ #include "compute_temperature.hh" #include "fft.hh" #include "material_point.hh" #include <cmath> ComputeTemperature::ComputeTemperature(Real dt) : dt(dt) {} /* -------------------------------------------------------------------------- */ void ComputeTemperature::setDeltaT(Real dt) { this->dt = dt; } /* -------------------------------------------------------------------------- */ void ComputeTemperature::compute(System& system) { UInt Nb = system.getNbParticles(); UInt dim = UInt(sqrt(Nb)); Matrix<complex> theta(dim); Matrix<complex> h_v(dim); Matrix<complex> theta_fourier(dim); Matrix<complex> h_v_fourier(dim); Matrix<std::complex<int>> wave_nb(dim); Matrix<complex> dthetadt(dim); // Initialize matrices for (auto&& par : index(theta)){ int i = std::get<0>(par); int j = std::get<1>(par); - UInt id = i*dim + j; + UInt id = i + j*dim; auto& ic = dynamic_cast<MaterialPoint&>(system.getParticle(id)); theta(i,j) = ic.getTemperature(); h_v(i,j) = ic.getHeatRate(); } theta_fourier = FFT::transform(theta); h_v_fourier = FFT::transform(h_v); wave_nb = FFT::computeFrequencies(dim); Real factor = 4.0*M_PI*M_PI/(L*L); for(int j=0; j<dim; j++){ for(int i=0; i<dim; i++){ dthetadt(i, j) = h_v_fourier(i, j) - k*theta_fourier(i, j)*factor*Real(norm(wave_nb(i, j))); dthetadt(i, j) = dthetadt(i, j)/(rho*C); } } dthetadt = FFT::itransform(dthetadt); for (auto&& par : index(dthetadt)){ int i = std::get<0>(par); int j = std::get<1>(par); - UInt id = i*dim + j; + UInt id = i + j*dim; auto& value = std::get<2>(par); dynamic_cast<MaterialPoint&>(system.getParticle(id)).getTemperature() += value.real()*this->dt; } } /* -------------------------------------------------------------------------- */ diff --git a/Homework3/src/compute_temperature.hh b/Homework3/src/compute_temperature.hh index 95373cc..4cb3eb3 100644 --- a/Homework3/src/compute_temperature.hh +++ b/Homework3/src/compute_temperature.hh @@ -1,33 +1,33 @@ #ifndef __COMPUTE_TEMPERATURE__HH__ #define __COMPUTE_TEMPERATURE__HH__ /* -------------------------------------------------------------------------- */ #include "compute.hh" //! Integrate heat equation class ComputeTemperature : public Compute { // Constructors/Destructors public: ComputeTemperature(Real timestep); // Methods public: //! Set time step void setDeltaT(Real dt); //! Evolve temperature void compute(System& system) override; private: Real dt; // variables are for now hard coded (-> to be modified) Real rho = 1.0; Real C = 1.0; Real k = 1.0; - Real L = 1.0; + Real L = 2.0; }; /* -------------------------------------------------------------------------- */ #endif //__COMPUTE_TEMPERATURE__HH__