diff --git a/Homework3/src/compute_temperature.cc b/Homework3/src/compute_temperature.cc index a61b81b..395c49c 100644 --- a/Homework3/src/compute_temperature.cc +++ b/Homework3/src/compute_temperature.cc @@ -1,84 +1,67 @@ #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); + Matrix<complex> dthetadt_fourier(dim); // Initialize matrices for (auto&& par : index(theta)){ int i = std::get<0>(par); int j = std::get<1>(par); UInt id = i + j*dim; auto& ic = dynamic_cast<MaterialPoint&>(system.getParticle(id)); theta(i,j) = complex(ic.getTemperature(),0); h_v(i,j) = complex(ic.getHeatRate() ,0); + + std::cout<<theta(i,j).real()<<std::endl; } - - for(int i=0; i<dim; i++) - { - for(int j=0; j<dim; j++) - { - std::cout<<theta(i,j).real()<<" "; - } - std::cout<<std::endl; - } - std::cout<<"========"<<std::endl; theta_fourier = FFT::transform(theta); h_v_fourier = FFT::transform(h_v); wave_nb = FFT::computeFrequencies(dim); - - for(int i=0; i<dim; i++) - { - for(int j=0; j<dim; j++) - { - std::cout<<theta_fourier(i,j).imag()<<" "; - } - std::cout<<std::endl; - } + Real factor = pow(2.0*M_PI/L,2); + for(int j=0; j<dim; j++){ for(int i=0; i<dim; i++){ - dthetadt(i, j) = -k*theta_fourier(i,j)*factor*Real(norm(wave_nb(i,j))) - +h_v_fourier(i,j); - dthetadt(i,j) = dthetadt(i,j)/(rho*C); + dthetadt_fourier(i, j) = (-k*theta_fourier(i,j)*factor*Real(norm(wave_nb(i,j))) + +h_v_fourier(i,j))/(rho*C); } } - dthetadt = FFT::itransform(dthetadt); + dthetadt = FFT::itransform(dthetadt_fourier); for (auto&& par : index(dthetadt)){ int i = std::get<0>(par); int j = std::get<1>(par); UInt id = i + j*dim; auto& value = std::get<2>(par); dynamic_cast<MaterialPoint&>(system.getParticle(id)).getTemperature() += value.real()*this->dt; } } /* -------------------------------------------------------------------------- */