diff --git a/homework3/starting_point/.idea/workspace.xml b/homework3/starting_point/.idea/workspace.xml
index 688e01b..ec5643a 100644
--- a/homework3/starting_point/.idea/workspace.xml
+++ b/homework3/starting_point/.idea/workspace.xml
@@ -1,461 +1,460 @@
-
-
+
-
-
+
+
-
+
-
+
-
+
-
-
+
+
-
-
+
+
-
+
-
+
-
+
true
DEFINITION_ORDER
1544171391000
1544171391000
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
+
-
-
-
-
-
-
+
+
-
+
-
-
+
+
-
+
-
-
+
+
+
+
+
+
-
+
-
-
+
+
-
+
-
+
-
+
-
-
+
+
\ No newline at end of file
diff --git a/homework3/starting_point/compute_temperature.cc b/homework3/starting_point/compute_temperature.cc
index 7d8ec26..da36890 100644
--- a/homework3/starting_point/compute_temperature.cc
+++ b/homework3/starting_point/compute_temperature.cc
@@ -1,77 +1,79 @@
#include "compute_temperature.hh"
#include "fft.hh"
#include "material_point.hh"
#include
/* -------------------------------------------------------------------------- */
ComputeTemperature::ComputeTemperature(Real dt) : dt(dt) {}
void ComputeTemperature::setDeltaT(Real dt) {
this->dt = dt;
}
/*Compute takes a list of particles creates matrices of T and h solve for T using spectral discr. in space and finite
differences in t*/
void ComputeTemperature::compute(System& system) {
UInt N = system.getNbParticles();
UInt Msize = sqrt(N);
Real K;
- Matrix Temp(N);
- Matrix hrate(N);
+ Matrix Temp(Msize);
+ Matrix hrate(Msize);
Real L = 2.;
- Real dx = L/N;
+ Real dx = L/Msize;
//Fill the matrix
for (auto & par : system) {
auto &mpoint = static_cast(par);
auto &T = mpoint.getTemperature();
auto &h = mpoint.getHeatRate();
auto &x = mpoint.getPosition();//vector x,y
- int i = int(x[0] / dx);
- int j = int(x[1] / dx);
+ int i = int((x[0]+L/2) / dx);
+ int j = int((x[1]+L/2) / dx);
Temp(i, j) = T;
hrate(i, j) = h;
}
+
//FFT forward
FFT matFFT;
auto Temp_sp = matFFT.transform(Temp);
auto hrate_sp = matFFT.transform(hrate);
auto freq = matFFT.computeFrequencies(Msize);
Matrix rhs_sp(Msize);
//compute right-hand-side in Spectral Space
auto q2 = abs2(freq);
Real fac = (K*4*M_PI*M_PI)/(L*L);
Temp_sp *= fac;
Temp_sp = Temp_sp*q2;
rhs_sp = hrate_sp -Temp_sp;
//FFT backward
auto rhs = matFFT.itransform(rhs_sp);
//update T
rhs *= dt;
Temp = Temp + rhs;
//update list of particles
for (auto & par : system) {
auto &mpoint = static_cast(par);
auto &x = mpoint.getPosition();//vector x,y
- int i = int(x[0] / dx);
- int j = int(x[1] / dx);
+ int i = int((x[0]+L/2) / dx);
+ int j = int((x[1]+L/2) / dx);
mpoint.getTemperature()= Temp(i,j).real();
}
+
}
/* -------------------------------------------------------------------------- */