StationaryHeatTransfer.py 1.24 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Feb  1 08:49:11 2019

@author: bleyerj
"""
from dolfin import *
import mfront_wrapper as mf
import numpy as np

length = 30e-3
width = 5.4e-3
mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)

16
V = FunctionSpace(mesh, "CG", 1)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
17 18 19 20 21 22 23 24 25 26 27 28 29
T = Function(V, name="Temperature")

def left(x, on_boundary):
    return near(x[1], 0) and on_boundary
def right(x, on_boundary):
    return near(x[0], length) and on_boundary

Tl = 300
Tr = 800

bc = [DirichletBC(V, Constant(Tl), left),
      DirichletBC(V, Constant(Tr), right)]

30
facets = MeshFunction("size_t", mesh, 1)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
31 32 33 34
ds = Measure("ds", subdomain_data=facets)



35 36
material = mf.MFrontNonlinearMaterial("materials/src/libBehaviour.so",
                                      "StationaryHeatTransfer",
Jeremy BLEYER's avatar
Jeremy BLEYER committed
37
                                      hypothesis="plane_strain")
38

39
problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
40
problem.bc = bc
41
# problem.register_gradient("TemperatureGradient", grad(T))
Jeremy BLEYER's avatar
Jeremy BLEYER committed
42

43
T.interpolate(Constant(Tl))
Jeremy BLEYER's avatar
Jeremy BLEYER committed
44 45 46

problem.solve(T.vector())

47 48 49 50 51
x = np.linspace(0, length)
import matplotlib.pyplot as plt
plt.plot(x, np.array([T(xi, width/2) for xi in x]))
#plt.figure()
#plot(project(flux.function[0], FunctionSpace(mesh, "DG",0)))