Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85081878
compute_temperature.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Sep 26, 16:08
Size
2 KB
Mime Type
text/x-c
Expires
Sat, Sep 28, 16:08 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21051152
Attached To
R9484 sp4e-homework-lars-bertil
compute_temperature.cc
View Options
#include "compute_temperature.hh"
#ifdef USE_FFTW
#include "fft.hh"
#endif
#include "material_point.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
void
ComputeTemperature
::
compute
(
System
&
system
)
{
#ifdef USE_FFTW
std
::
cout
<<
density
<<
std
::
endl
;
UInt
size
=
system
.
getNbParticles
();
UInt
N
=
std
::
sqrt
(
size
);
Matrix
<
complex
>
temperature
(
N
);
Matrix
<
complex
>
heat_source
(
N
);
std
::
transform
(
system
.
begin
(),
system
.
end
(),
temperature
.
begin
(),
[](
auto
&
part
)
{
return
static_cast
<
MaterialPoint
&>
(
part
).
getTemperature
();
});
std
::
transform
(
system
.
begin
(),
system
.
end
(),
heat_source
.
begin
(),
[](
auto
&
part
)
{
return
static_cast
<
MaterialPoint
&>
(
part
).
getHeatSource
();
});
auto
freqs
=
FFT
::
computeFrequencies
(
N
);
auto
h_hat
=
FFT
::
transform
(
heat_source
);
auto
T_hat
=
FFT
::
transform
(
temperature
);
auto
T_rate_hat
=
Matrix
<
complex
>
(
N
);
const
auto
factor
=
delta_t
*
density
/
capacity
;
const
auto
factor2
=
1.
/
density
/
capacity
;
auto
implicit_integrator
=
[
&
](
auto
&&
entry
)
{
int
i
=
std
::
get
<
0
>
(
entry
);
int
j
=
std
::
get
<
1
>
(
entry
);
auto
&
T_rate
=
std
::
get
<
2
>
(
entry
);
auto
f
=
freqs
(
i
,
j
);
Real
k2
=
std
::
abs
(
complex
(
f
.
real
(),
f
.
imag
()))
*
2.
*
M_PI
/
L
;
k2
*=
k2
;
const
auto
factor2
=
1.
+
factor
*
conductivity
*
k2
;
T_rate_hat
(
i
,
j
)
=
1.
/
factor2
*
(
factor
*
h_hat
(
i
,
j
)
+
T_hat
(
i
,
j
))
-
T_hat
(
i
,
j
);
};
auto
explicit_integrator
=
[
&
](
auto
&&
entry
)
{
int
i
=
std
::
get
<
0
>
(
entry
);
int
j
=
std
::
get
<
1
>
(
entry
);
auto
&
T_rate
=
std
::
get
<
2
>
(
entry
);
auto
f
=
freqs
(
i
,
j
);
Real
k2
=
std
::
abs
(
complex
(
f
.
real
(),
f
.
imag
()))
*
2.
*
M_PI
/
L
;
k2
*=
k2
;
T_rate_hat
(
i
,
j
)
=
1.
/
factor2
*
(
h_hat
(
i
,
j
)
+
T_hat
(
i
,
j
)
*
conductivity
*
k2
);
};
if
(
implicit
)
for
(
auto
entry
:
index
(
T_rate_hat
))
implicit_integrator
(
entry
);
else
for
(
auto
entry
:
index
(
T_rate_hat
))
explicit_integrator
(
entry
);
auto
temperature_rate
=
FFT
::
itransform
(
T_rate_hat
);
auto
it_temp_rate
=
temperature_rate
.
begin
();
for
(
auto
&
part
:
system
)
{
auto
&
p
=
static_cast
<
MaterialPoint
&>
(
part
);
p
.
getTemperatureRate
()
=
(
*
it_temp_rate
).
real
();
p
.
getTemperature
()
+=
delta_t
*
p
.
getTemperatureRate
();
++
it_temp_rate
;
}
#endif
}
/* -------------------------------------------------------------------------- */
Event Timeline
Log In to Comment