Heat-Equation/Readme.md

```python
#1 Loading functions  and modules
from fenics import *
import matplotlib.pyplot as plt
T = 0.1
num_steps = 20
dt = T/num_steps
rho = 7500
Cp = 500
k = 50
alpha = k/(rho*Cp)
```


```python
#2 Create mesh and define function space
nx = 0.008
ny = 0.003
mesh = RectangleMesh(Point(0,0),Point(nx,ny),30, 30,'left')
V = FunctionSpace(mesh, 'Lagrange', 1) #Lagrange are triangular elements
plot(mesh)
plt.show()

```


![png](output_1_0.png)



```python
# Boundary conditions
u0 = Constant(100)
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V,u0, boundary)
```


```python
u_n = project(1, V)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
F = u*v*dx + alpha*dt*dot(grad(u), grad(v))*dx-u_n*v*dx
a, L = lhs(F), rhs(F)
```


```python
vtkfile = File('solution/solution.pvd')
u = Function(V)
t = 0
for n in range(num_steps):
    t += dt
    #u0.t = t
    solve(a == L, u, bc)
    #c = plot(u,)
    #plt.colorbar(c)
    #plt.show()
    ####
    vtkfile << (u, t)
    u_n.assign(u)
```


```python
1E-13
```




    1e-13




```python
1e2
```




    100.0




```python
1E-13+1e2
```




    100.0000000000001



# Boundary conditions


```python
#1 Loading functions  and modules
from fenics import *
import matplotlib.pyplot as plt
T = 4
num_steps = 200
dt = T/num_steps
rho = 7500
Cp = 500
k = 50
alpha = k/(rho*Cp)
```


```python
#2 Create mesh and define function space
nx = 0.008
ny = 0.003
mesh = RectangleMesh(Point(0,0),Point(nx,ny),8, 8,'left')
V = FunctionSpace(mesh, 'Lagrange', 1) #Lagrange are triangular elements
plot(mesh)
plt.show()


```


![png](output_10_0.png)



```python
#Boundary Conditions
tol = 1E-14

def BC1(x, on_boundary):
    return on_boundary and abs(x[0]-nx) < tol

def BC2(x, on_boundary):
    return on_boundary and abs(x[0]-0) < tol

bc1=DirichletBC(V,Constant(25),BC1)
bc2=DirichletBC(V,Constant(800),BC2)
bc=(bc1,bc2)
```


```python
u_n = project(25, V)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
F = u*v*dx + alpha*dt*dot(grad(u), grad(v))*dx-u_n*v*dx
a, L = lhs(F), rhs(F)
```


```python
vtkfile = File('solution/solution2.pvd')
u = Function(V)
t = 0
for n in range(num_steps):
    t += dt
    #u0.t = t
    solve(a == L, u, bc)
    vtkfile << (u, t)
    #c = plot(u,)
    #plt.colorbar(c)
    #plt.show()
    u_n.assign(u)
```


```python

```