Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88767405
superposition.py
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
Sun, Oct 20, 14:20
Size
4 KB
Mime Type
text/x-python
Expires
Tue, Oct 22, 14:20 (2 d)
Engine
blob
Format
Raw Data
Handle
21808628
Attached To
R11910 Additive Manufacturing Work
superposition.py
View Options
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Dec 22 16:08:08 2021
@author: kubilay
"""
from
fenics
import
*
import
os
import
numpy
as
np
from
helpers
import
*
import
time
as
timer
import
multiprocessing
as
mp
#define simulation parameters
velocity
=
0.8
*
np
.
array
([
1
,
0
,
0
])
velocity_mag
=
np
.
linalg
.
norm
(
velocity
)
source_loc
=
np
.
array
([
-
1.95
,
-
1.95
,
4
])
alpha
=
3
time
=
1
V
=
velocity
/
(
2
*
alpha
)
N
=
20
num_steps
=
40
dt
=
time
/
num_steps
tol
=
1e-14
#read the mesh
mesh
=
Mesh
(
'../Part_geometry/mesh/layer_004.xml'
)
boundary_mesh
=
BoundaryMesh
(
mesh
,
'exterior'
)
#determine the limits of the domain
x_max
=
np
.
max
(
mesh
.
coordinates
()[:,
0
])
x_min
=
np
.
min
(
mesh
.
coordinates
()[:,
0
])
y_max
=
np
.
max
(
mesh
.
coordinates
()[:,
1
])
y_min
=
np
.
min
(
mesh
.
coordinates
()[:,
1
])
z_max
=
np
.
max
(
mesh
.
coordinates
()[:,
2
])
z_min
=
np
.
min
(
mesh
.
coordinates
()[:,
2
])
#define function space for v
V
=
FunctionSpace
(
mesh
,
'P'
,
1
)
#define markers for the boundary faces
boundary_markers
=
MeshFunction
(
'size_t'
,
mesh
,
mesh
.
topology
()
.
dim
()
-
1
,
4
)
#define boundary subclasses for bottom, top and side surfaces, mark the boundaries
#define boundary subclasses for bottom, top and side surfaces, mark the boundaries
class
BoundaryBottom
(
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
on_boundary
and
near
(
x
[
2
],
z_min
,
tol
)
b_bottom
=
BoundaryBottom
()
b_bottom
.
mark
(
boundary_markers
,
0
)
class
BoundaryTop
(
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
on_boundary
and
near
(
x
[
2
],
z_max
,
tol
)
b_top
=
BoundaryTop
()
b_top
.
mark
(
boundary_markers
,
1
)
class
BoundarySides
(
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
on_boundary
and
(
near
(
x
[
0
],
x_min
,
tol
)
or
near
(
x
[
0
],
x_max
,
tol
)
or
near
(
x
[
1
],
y_min
,
tol
)
or
near
(
x
[
1
],
y_max
,
tol
))
b_sides
=
BoundarySides
()
b_sides
.
mark
(
boundary_markers
,
2
)
#initial source location
class
BoundarySource
(
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
on_boundary
and
(
near
(
x
[
2
],
z_max
,
tol
)
and
((
x
[
0
]
-
1.8
)
**
2
+
x
[
1
]
**
2
<
0.14
))
b_source
=
BoundarySource
()
b_source
.
mark
(
boundary_markers
,
3
)
#redefine d in therm of the boundary markers
ds
=
Measure
(
'ds'
,
domain
=
mesh
,
subdomain_data
=
boundary_markers
)
#define initial condition
T_d
=
Constant
(
0
)
#define Neumann boundary conditions
g_top
=
Constant
(
0
)
g_top
=
Constant
(
0.0
)
g_sides
=
Constant
(
0
)
source
=
Constant
(
0
)
#list of boundary conditions for convenience
boundary_conditions
=
{
0
:
{
'Dirichlet'
:
T_d
},
1
:
{
'Neumann'
:
g_top
},
2
:
{
'Neumann'
:
g_sides
},
3
:
{
'Neumann'
:
source
}}
#project the initial condition
T_0
=
project
(
Constant
(
0.0
),
V
)
#define test and trial function and the source term
T
=
TrialFunction
(
V
)
v
=
TestFunction
(
V
)
f
=
Constant
(
0
)
#gather all Dirichlet boundary conditions in one list
bcs
=
[]
for
i
in
boundary_conditions
:
if
'Dirichlet'
in
boundary_conditions
[
i
]:
bc
=
DirichletBC
(
V
,
boundary_conditions
[
i
][
'Dirichlet'
],
boundary_markers
,
i
)
bcs
.
append
(
bc
)
#gather all Neumann boundary conditions in weak form
integrals_N
=
[]
for
i
in
boundary_conditions
:
if
'Neumann'
in
boundary_conditions
[
i
]:
if
boundary_conditions
[
i
][
'Neumann'
]
!=
0
:
g
=
boundary_conditions
[
i
][
'Neumann'
]
integrals_N
.
append
(
g
*
v
*
ds
(
i
))
#write the weak form and seperate into right- and left-hand sides
F
=
T
*
v
*
dx
+
dt
*
dot
(
grad
(
T
),
grad
(
v
))
*
dx
-
T_0
*
v
*
dx
-
dt
*
f
*
v
*
dx
+
dt
*
sum
(
integrals_N
)
a
,
L
=
lhs
(
F
),
rhs
(
F
)
#define the function and time
T
=
Function
(
V
)
t
=
0
#create directory to save files
try
:
os
.
mkdir
(
'results'
)
except
FileExistsError
:
for
file
in
os
.
scandir
(
'results'
):
os
.
remove
(
file
.
path
)
#define vtk file to save the data at each iteration
vtkfile
=
File
(
'results/analytic_temperature.pvd'
)
def
run_analytic_analysis
(
t
):
analytical_temperature
=
MyExpression
(
time
=
t
,
velocity
=
velocity
,
source_loc
=
source_loc
+
velocity
*
t
,
alpha
=
alpha
)
T
=
interpolate
(
analytical_temperature
,
V
)
print
(
np
.
max
(
T
.
vector
()[:]))
T
.
rename
(
"Temperature"
,
"Analytic Temperature"
)
vtkfile
<<
(
T
,
t
)
start
=
timer
.
time
()
#T_final = Function(V)
#solve the weak form for each timestep
for
n
in
range
(
num_steps
):
#increment time and recalculate time dependent functions
t
+=
dt
T_d
.
t
=
t
g_top
.
t
=
t
g_sides
.
t
=
t
run_analytic_analysis
(
t
)
print
(
'it took:'
,
timer
.
time
()
-
start
)
Event Timeline
Log In to Comment