Commit 4e88d175 authored by Remi Lacombe's avatar Remi Lacombe
Browse files

initial version

parents
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np
# Mesh + Function Space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'Lagrange', 1) # 'DP' for discontinuous Galerkin
# Pour un mesh sur le cercle unite
''' domain = Circle(Point(0, 0), 1)
mesh = generate_mesh(domain, 64) '''
# Boundary Condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2) # degree donne le degre d'interpolation de cette expression dans chaque element
# x[0] est la coordonnee x et x[1] est la coordonnee y
def boundary(x, on_boundary): # on_boundary is a boolean
return on_boundary
# on_boundary est donne par FEniCS pour tout point x du domaine. On peut laisser comme ca ou bien definir une autre frontiere explicitement
# Dans ce cas attention a bien utiliser une tolerance pour eviter les erreurs d'arrondis ! La fonction near() peut etre utilisee pour cela
bc = DirichletBC(V, u_D, boundary) # boundary is the subdomain here
# Variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx # dot est le produit scalaire. Utiliser inner pour des matrices tho. Definit le PDE a resoudre
L = f*v*dx
# Pour interpoler une expression (fonction continue) en une fonction discrete:
''' ff = interpolate(f, V) '''
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Plot
#plot(u, title='Finite element solution')
#plot(mesh, title='Finite element mesh')
# Save solution VTK format
vtkfile = File('poisson/solution.pvd')
vtkfile << u
# Error in L2 norm
error_L2 = errornorm(u_D, u, 'L2')
# Compute maximum error at vertice
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh) # Renvoit un tableau avec le valeurs de u
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))
print('Erreur L = ' , error_L2)
print('Erreur maximale = ', error_max)
# Convertit dans un tableau les valeurs de u sur les noeuds du maillage. Attention l'ordre est different e celui de la fonction precedente
nodal_values_u = u.vector()
array_u = nodal_values_u.array()
interactive() # Necessaire pour maintenir le plot.
# Aide pour le trace de fonctions
''' tol = 0.001
y = np.linspace(-1 + tol, 1 - tol, 101)
points = [(0, y_) for y_ in y] #points d'abscisse 0
w_line = np.array([w(point) for point in points])
plt.plot(y, w_line, 'k', linewidth=2) # use "b--" instead of "k" pour pointilles
plt.grid(True)
plt.xlabel('$y$')
plt.legend(['Deflection', 'Load'], loc = 'upper left')
plt.savefig('poisson/curves.pdf')
plt.savefig('poisson/curves.png') '''
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment