demo multi-material

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

# -*- coding: utf-8 -*-
"""
Problema de 2a.-ordem unidimensional
-(a * u,x),x - q = 0, x in ]0,1[
Condição de Dirichlet
u(0) = 1
u(1) = 1
"""
import matplotlib.pyplot as plt
from dolfin import *
# definição da malha para problema unidimensional
mesh = IntervalMesh(50, 0, 1)
# Definição do espaço de função na malha utilizando polinômios de Lagrange do 1o. grau
V = FunctionSpace(mesh, "Lagrange", 1)
# Condição de Dirichlet à esquerda
def LeftDirichletBoundary( x, on_boundary):
 # Função que verifica se ponto x pertence ao LeftDirichletBoundary com tolerânica de +/- tol
 # retorna veradeiro ou falso
 tol = 1E-14
 # tolerância necessária ao lidar com representação de número em ponto flutuante
 return on_boundary and abs(x[0]) < tol
# Condição de Dirichlet à direita
def RightDirichletBoundary( x, on_boundary):
 # Função que verifica se ponto x pertence ao LeftDirichletBoundary com tolerânica de +/- tol
 # retorna veradeiro ou falso
 tol = 1E-14
 # tolerância necessária ao lidar com representação de número em ponto flutuante
 return on_boundary and abs(x[0] - 1) < tol
# Definição do valor no contorno esquerdo
u_L = Constant(0) 
# Definição do valor no contorno direito
u_R = Constant(1)
# Imposição das condições de contorno de Dirichlet à esquerda e à direita
bcs = [DirichletBC(V, u_L, LeftDirichletBoundary), DirichletBC(V, u_R, RightDirichletBoundary)]
# Definição da forma fraca
u = TrialFunction(V) # função de aproximação
v = TestFunction(V) # função peso
q = Constant(1) # variável independente
#a = Expression('x[0]<=0.5+tol ? x[0] : 2*(1-x[0])',degree=1,tol=1e-14) # 2 multi-material
a = Expression('x[0]<=0.5+tol ? x[0] : x[0]<=0.75+tol ? .010 : 2',degree=1,tol=1e-14) # 3 multi-material
B = a*inner(grad(u), grad(v))*dx # Operador bilinear
L = q*v*dx # operador linear
# Resolve o problema
u = Function(V)
solve( B == L, u, bcs)
# Plota a solução
plot(u,label='2 elementos')
plt.legend()
plt.title('Campo de temperatura')
# Determina u em x=0.5
print(u(0.5))