<jupyter_start><jupyter_text>Basic Poisson equationThis example is implemented in the Python file poisson_basic.py and it illustrates how to:- Solve a linear partial differential equation;- Create a mesh for a disk domain;- Define a `FunctionSpace`;- Create and apply Dirichlet boundary conditions;- Plot the solution;- Use `errornorm` to compare results with the exact solution;- Evaluate the value of the function at an arbitrary location.  Equation and problem definition The Poisson equation on a unit disk ($R=1$) with  Dirichlet boundary conditions:        $\begin{array}{r c c l}    \nabla^2 u &=& -1 & \text{in } \Omega,\\    u &=& 0 & \text{on } \Gamma.\\    \end{array}$    &emsp;        The exact solution in Cartesian coordinates is$$u({\bf x})\equiv u(x,y)=\frac{1}{4} \left(1-x^2-y^2\right).$$  Weak formulation of the problem In FEniCS, we are solving the weak form of the PDE$$\int_\Omega d{\bf x} \ (\nabla^2 u + 1) \ v = 0,$$where $v$ is an arbitrary *test function* that evaluates to $0$ on the boundary points with prescribed values (Dirichlet boundary conditions). Using integration by parts, the above equation can be rewritten as $$\int_\Omega d{\bf x} \ (-\nabla u \cdot \nabla v + v) = 0.$$  Implementation In order to use the FEniCS package we need to import relevant libraries:<jupyter_code>from __future__ import print_function
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>Define a circular domain and generate a triangular mesh. The cell size is approximately equal to the diameter of the domain divided by the resolution ($N$).<jupyter_code># Create mesh
R = 1. # radius
N = 20 # mesh resolution
domain = Circle(Point(0., 0.), R)
mesh = generate_mesh(domain, N)
plot(mesh,linewidth=0.3)
plt.show()<jupyter_output><empty_output><jupyter_text>Next, we define the `FunctionSpace` that lives on the mesh.<jupyter_code>#define function space
degreeElements = 1
FS = FunctionSpace(mesh, 'Lagrange', degreeElements)<jupyter_output><empty_output><jupyter_text>In this example, we use linear *Lagrange* elements. Higher order elements can be used by changing the value of the `degreeElements`.  The list of finite element types that are supported in FEniCS is available here: click this link. To impose Dirichlet boundary conditions on the `FunctionSpace`, we have to define a function that returns a boolean value: `True` if the given point $\bf x$ lies on the Dirichlet boundary and `False`otherwise. The argument `on_boundary` is `True` if $\bf x$ is on the physical boundaryof the mesh. Note that the boundary function will be called for every discrete pointin the mesh, which means that we may define boundaries also inside the physical domain.<jupyter_code>#impose Dirichlet boundary conditions
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(FS, Constant(0.), boundary)<jupyter_output><empty_output><jupyter_text>Define function $u(\bf x)$ and test function $v(\bf x)$.<jupyter_code>#define function u and test function v
u = Function(FS)
v = TestFunction(FS)<jupyter_output><empty_output><jupyter_text>The weak formulation of the problem$$\int_\Omega d{\bf x} \ (-\nabla u \cdot \nabla v + v) = 0,$$where `dx` describes the integration over the domain $\Omega$<jupyter_code># weak formulation of the problem
Res = -dot(grad(u), grad(v))*dx + v*dx<jupyter_output><empty_output><jupyter_text>Solve the problem using the Dirichlet boundary conditions `bc` and store the solution in the function `u`.<jupyter_code># solve the problem
solve(Res == 0, u, bc)<jupyter_output><empty_output><jupyter_text>Note that the `solve` function can be used to solve both linear and non-linear problems. The default solver is using the Newton-Raphson method. Consult with the  FEniCS tutorial for other solver options. The solution $u({\bf x})$ can be plotted with the `plot` option, which uses the `matplotlib` package.<jupyter_code># plot solution
c = plot(u,mode='color',title='$u$')
plt.colorbar(c)
plot(mesh,linewidth=0.3)
plt.show()<jupyter_output><empty_output><jupyter_text>The accuracy of the obtained solution $u({\bf x})$ can be compared against the exact solution$u({\bf x})=\frac{1}{4} \left(1-x^2-y^2\right)$ by calculating the $L^2$ norm of the difference:<jupyter_code># exact solution
uExact=Expression('(1-x[0]*x[0]-x[1]*x[1])/4',degree=2)

# Compute error (L2 norm)
error_L2 = errornorm(uExact, u, 'L2')
print("Error = ",error_L2)<jupyter_output>Error =  0.0010054862645973131<jupyter_text>The `Expression` is using the C++ syntax to define functions, where `x[0]` and `x[1]` refer to the $x$ and $y$ coordinates, respectively, and the `degree` argument specifies the order of polynomial functions used to interpolate the expression on mesh elements.  Function can be evaluated at an arbitrary point on the physical domain:<jupyter_code># evaluate function at several points
print("u(0,0) = ",u(0,0))
print("u(0.5,0.5) = ",u(0.5,0.5))<jupyter_output>u(0,0) =  0.2497371456650515
u(0.5,0.5) =  0.12449260452913821<jupyter_text>Note that function cannot be evaluated outside the physical domain<jupyter_code>print("u(2,0) = ",u(2,0))<jupyter_output><empty_output><jupyter_text>Be careful with points $x$ on the boundary $\Gamma$, because they may lie outside the discretized domain<jupyter_code>print("u(1,0) = ",u(1,0))
print("u(0,1) = ",u(0,1))<jupyter_output>u(1,0) =  0.0<jupyter_text>Plot the value of function $u(x,0)$ along the $x$ coordinate axis using the `matplotlib` functions.<jupyter_code># plot solution
tol=1e-2;
x=np.linspace(-1+tol,1-tol,40)
points = [(x_, 0) for x_ in x]
u_line = np.array([u(point) for point in points])
plt.plot(x,u_line,'k-')
plt.plot(x,(1-x*x)/4,'r--')
plt.xlabel('$x$')
plt.ylabel('$u(x,0)$')
plt.legend(['FEM','exact'])
plt.show()<jupyter_output><empty_output><jupyter_text>Complete code<jupyter_code>from __future__ import print_function
from fenics import *
from mshr import *
import matplotlib.pyplot as plt
import numpy as np

# Create mesh
R = 1. # radius
N = 20 # mesh resolution
domain = Circle(Point(0., 0.), R)
mesh = generate_mesh(domain, N)
plot(mesh,linewidth=0.3)
plt.show()

#define function space
degreeElements = 1
FS = FunctionSpace(mesh, 'Lagrange', degreeElements)

#impose Dirichlet boundary conditions
def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(FS, Constant(0.), boundary)

#define function u and test function v
u = Function(FS)
v = TestFunction(FS)

# weak formulation of the problem
Res = -dot(grad(u), grad(v))*dx + v*dx

# solve the problem
solve(Res == 0, u, bc)

# plot solution
c = plot(u,mode='color',title='$u$')
plt.colorbar(c)
plot(mesh,linewidth=0.3)
plt.show()

# exact solution
uExact=Expression('(1-x[0]*x[0]-x[1]*x[1])/4',degree=2)

# Compute error (L2 norm)
error_L2 = errornorm(uExact, u, 'L2')
print("Error = ",error_L2)

# evaluate function at several points
print("u(0,0) = ",u(0,0))
print("u(0.5,0.5) = ",u(0.5,0.5))
#print("u(2,0) = ",u(2,0)) #point outside domain
print("u(1,0) = ",u(1,0))
# print("u(0,1) = ",u(0,1)) #point outside the discretized domain


# plot solution
tol=1e-2;
x=np.linspace(-1+tol,1-tol,40)
points = [(x_, 0) for x_ in x]
u_line = np.array([u(point) for point in points])
plt.plot(x,u_line,'k-')
plt.plot(x,(1-x*x)/4,'r--')
plt.xlabel('$x$')
plt.ylabel('$u(x,0)$')
plt.legend(['FEM','exact'])
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1372284442045608, 'fraction_numerical': 0.027391715018216165, 'mean_word_length': 5.576752440106477, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 54, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '557833', 'license': 'MIT', 'n_tokens_mistral': 2644, 'n_tokens_neox': 2497, 'n_words': 825, 'path': 'poisson/poisson_basic.ipynb', 'repo_name': 'gletreut/FEniCS_tutorial', 'size': '7411'}