Page MenuHomec4science
Files F78055942

tutorial.py
No OneTemporary
Actions

File Metadata

Created
Sun, Aug 18, 03:33

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

c4science · Help