Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F83496226
tutorial.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
Tue, Sep 17, 11:22
Size
4 KB
Mime Type
text/x-python
Expires
Thu, Sep 19, 11:22 (2 d)
Engine
blob
Format
Raw Data
Handle
20810994
Attached To
R11910 Additive Manufacturing Work
tutorial.py
View Options
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Dec 3 12:35:17 2021
@author: ekinkubilay
"""
#from dolfin import *
#from dolfin_utils import meshconvert
from
fenics
import
*
import
os
import
numpy
as
np
#define simulation parameters
alpha
=
3
beta
=
1.2
time
=
10
num_steps
=
40
dt
=
time
/
num_steps
tol
=
1e-8
velocity
=
0.36
#read the mesh
mesh
=
Mesh
(
'../Part_geometry/mesh/layer_004.xml'
)
#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
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
=
Expression
(
'10*(-4+(x[1]*x[1]-2*t)+alpha*x[0]*x[0])'
,
degree
=
2
,
alpha
=
alpha
,
t
=
0
)
g_top
=
Constant
(
0.0
)
g_sides
=
Constant
(
0
)
source
=
Constant
(
-
500
)
#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/heat_equation_tutorial.pvd'
)
#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
#solve the weak form
solve
(
a
==
L
,
T
,
bcs
)
#save the temperature data the the .vtk file
vtkfile
<<
(
T
,
t
)
#assign T as the initial condition for the next iteration
T_0
.
assign
(
T
)
#remove Source Boundary from markers and RHS
b_source
.
mark
(
boundary_markers
,
4
)
L
-=
dt
*
g
*
v
*
ds
(
3
)
#redefine Source Boundary with new position and mark
class
BoundarySource
(
SubDomain
):
def
inside
(
self
,
x
,
on_boundary
):
return
on_boundary
and
(
near
(
x
[
2
],
z_max
,
tol
)
and
((
x
[
0
]
-
1.8
+
velocity
*
t
)
**
2
+
x
[
1
]
**
2
<
0.14
))
b_source
=
BoundarySource
()
b_source
.
mark
(
boundary_markers
,
3
)
#redefine d in therm of the new boundaries
ds
=
Measure
(
'ds'
,
domain
=
mesh
,
subdomain_data
=
boundary_markers
)
#add new Neumann boundary (Source Boundary) to RHS
L
+=
dt
*
g
*
v
*
ds
(
3
)
Event Timeline
Log In to Comment