<jupyter_start><jupyter_text>Quick overview of the finite element method -->Example on a complicated domain for solving PDEs.The finite element method is a rich and versatile approach to constructcomputational schemes to solve any partial differential equation onany domain in any dimension. The method may at first glance appearcumbersome and even unnatural as it relies on variational formulationsand polynomial spaces.Let us start by outlining the concepts briefly.Consider the following PDE in 2D: $$-\nabla^2 u = -u_{xx} - u_{yy} = f,$$ equipped with suitable boundary conditions.A finite difference scheme to solve the current PDEwould in the simplest case be described by the stencil $$\begin{equation}\label{overview:2d:fdm0} \tag{1}-\frac{u_{i-1,j} - 2 u_{i,j} + u_{i+1,j}}{h^2}-\frac{u_{i,j-1} - 2 u_{i,j} + u_{i,j+1}}{h^2} = f_{i}\end{equation}$$ or reordered to the more recognizable $$\begin{equation}\label{overview:2d:fdm} \tag{2}\frac{-u_{i-1,j} -u_{i,j-1} + 4 u_{i,j} - u_{i+1,j} -u_{i,j+1}}{h^2}  = f_{i}{\thinspace .}\end{equation}$$ On a structured mesh, the stencil appears natural andis convenient to implement.However, for a unstructured, "complicated" domainas shown in [Figure](overview:meshex),we would need to be careful when placingpoints and evaluating stencils and functions.   In particular, it will be difficult to evaluate the stencil near the dolphin in [Figure](overview:meshex) because some points will be on the inside and some outside on the outside of the dolphin. Both accuracy and efficiencymay easily be sacrificed by a reckless implementation.In general, a domain like the one represented in [Figure](overview:meshex) will be represented by a triangulation. Thefinite element method (and the finite volume method which often is aspecial case of the finite element method) is a methodology forcreating stencils in a structured mannerthat adapt to the underlying triangulation.The triangulation in [Figure](overview:meshex) is a mesh thatconsists of cells that are connected and defined in terms ofvertices. The fundamental idea of the finite element method isto construct a procedure to compute a stencil on a general element andthen apply this procedure to each element of the mesh. Letus therefore denote the mesh as $\Omega$ while $\Omega_e$ is the domainof a generic element such that $\Omega=\cup_e \Omega_e$.This is exactly the point where the challenges of the finite elementmethod start and where we need some new concepts.  The basic questionis: How should we create a stencil for ageneral element and a general PDE that has the maximal accuracy andminimal computational complexity at the current triangulation?  Thetwo basic building blocks of the finite element method are1. the solution is represented in terms of a polynomial expression on the   given general element, and2. a variational formulation of the PDE   where element-wise integration enables the PDE to be transformed to a   stencil.Step 1 is, as will be explained later, conveniently representedboth implementation-wise and mathematically as a solution $$\begin{equation}\label{overview:u:fem} \tag{3}u = \sum_{i=0}^N c_i {\psi}_i(x,y),\end{equation}$$ where $\{c_i\}$ are the coefficients to be determined(often called the degrees of freedom)and ${\psi}_i(x,y)$ are prescribed polynomials.The basis functions ${\psi}_i(x,y)$ used to express the solutionis often called the trial functions.  The next step is the variational formulation. This stepmay seem like a magic trick or a cumbersomemathematical exercise at first glance.We take the PDE and multiply by a function $v$ (usually called the test function)and integrate over an element $\Omega_e$ and obtain the expression $$\begin{equation}\label{overview:poisson} \tag{4}\int_{\Omega_e} -\nabla^2 u \, v {\, \mathrm{d}x} =  \int_{\Omega_e} f \, v {\, \mathrm{d}x}\end{equation}$$ A perfectly natural question at this point is: Why multiplywith a test function $v$? The simple answer is thatthere are $N+1$ unknowns that need to be determined in $u$in ([3](overview:u:fem))and for this we need $N+1$ equations. The equations areobtained by using  $N+1$ different test functions which when usedin  ([5](overview:fem:a))give rise to $N+1$ linearly independent equations.While ([4](overview:poisson)) is a variational formulation ofour PDE problem, it is not the most common form.It is common to re-write $$\begin{equation}\label{overview:fem:a} \tag{5}\int_{\Omega_e} -\nabla^2 u \, v {\, \mathrm{d}x}\end{equation}$$ to weaken the requirement of the polynomial space used for thetrial functions (that here needs to be twice differentiable)and write this term in its corresponding weak form.Thatis, the term is rewritten in terms of first-derivatives only (ofboth the trial and the test function)  with the aid of Gauss-Green's lemma: $$\begin{equation}\label{overview:fem:a:weak} \tag{6}\int_{\Omega_e} -\nabla^2 u \, v {\, \mathrm{d}x}  =\int_{\Omega_e} \nabla u \cdot \nabla v {\, \mathrm{d}x}   - \int_{\partial \Omega_e} \frac{\partial u}{\partial n} \,  v \,  dS\end{equation}$$ The reasons behind this alternative formulation are rather mathematical and willnot be a major subject of this book as they are well described elsewhere.In fact, a precise explanation would need tools from functional analysis.With the above rewrite and assuming now that the boundary term vanishes due toboundary conditions (why this is possible will be dealt with in detaillater) thestencil, corresponding to ([2](overview:2d:fdm)),  is represented by $$\int_{\Omega_e} \nabla u \cdot \nabla v {\, \mathrm{d}x}$$ where $u$ is called the *trial function*, $v$ is called a *test function*,and $\Omega$ is an element ofa triangulated mesh. The idea of software like FEniCS is that thispiece of mathematics can be directly expressed in terms of Python code as<jupyter_code># DO NOT RUN THIS CELL, THIS IS A DEMONSTRATION
mesh = Mesh("some_file")
V = FunctionSpace(mesh, "some polynomial")
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx<jupyter_output><empty_output><jupyter_text>The methodology and code in this example is not tied to a particularequation, except the formula for `a`, holding the derivatives of oursample PDE, but any other PDE terms could be expressed via `u`, `v`,`grad`, and other symbolic operators in this line of code.  In fact,finite element packages like FEniCS are typically structured asgeneral toolboxes that can be adapted to any PDE as soon as thederivation of variational formulations is mastered.  The main obstaclehere for a novice FEM user is then to understand the concept of trialfunctions and test functions realized in terms of polynomial spaces.Hence, a finite element formulation (or a weak formulation) ofthe Poisson problem that works on any mesh $\Omega$ can be writtenin terms of solving the problem: $$\int_\Omega\nabla u\cdot\nabla vd{\, \mathrm{d}x} = \int_\Omega fv{\, \mathrm{d}x}{\thinspace .}$$ By varying the trial and test spaces we obtain different stencils,some of which will be identical to finite difference schemes onparticular meshes. We will now show a complete FEniCS program toillustrate how a typical finite element code may be structured<jupyter_code># DO NOT RUN THIS CELL, THIS IS A DEMONSTRATION
mesh = Mesh("some_file")
V = FunctionSpace(mesh, "some polynomial")
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

bc = DirichletBC(V, "some_function", "some_domain")
solution = Function(V)  # unknown FEM function
solve(a == L, solution, bc)
plot(solution)<jupyter_output><empty_output><jupyter_text>While the finite element method is versatile and may be adapted to anyPDE on any domain in any dimension, the different methods that arederived by using different trial and test functions may varysignificantly in terms of accuracy and efficiency. In fact, a bad choice of polynomial space may in some cases lead to acompletely wrong result. This is particularly the case for complicatedPDEs. For this reason, it is dangerous to regard the method as a blackbox and not do proper verification of the method for a  particularapplication.In our view, thereare three important tests that should be frequently employedduring verification:1. reducing the model problem to 1D and carefully check the calculations involved in the variational formulation on a small 1D mesh2. perform the calculation involved on one general or random element3. test whether convergences is obtained and to what order the method converge by refining the meshThe two first tasks here should ideally be performed by independent calculationsoutside the framework used for the simulations. In our view `sympy` is aconvenient tool that can be used to assist hand calculations.So far, we have outlined how the finite element method handles derivativesin a PDE, but we also had a right-hand side function $f$. This term is multipliedby the test function $v$ as well, such that the entire Poisson equationis transformed to $$\int_\Omega\nabla u\cdot\nabla vd{\, \mathrm{d}x} = \int_\Omega fv{\, \mathrm{d}x}{\thinspace .}$$ This statement is assumed valid for all test functions $v$ in somefunction space $V$ of polynomials. The right-hand side expression iscoded in FEniCS as<jupyter_code># DO NOT RUN THIS CELL. RUN THE FULL PROGRAM (IN THE END) INSTEAD.
L = f*v*dx<jupyter_output><empty_output><jupyter_text>and the problem is then solved by the statements<jupyter_code># DO NOT RUN THIS CELL. RUN THE FULL PROGRAM (IN THE END) INSTEAD.
u = Function(V)  # unknown FEM function
solve(a == L, u, bc)<jupyter_output><empty_output><jupyter_text>where `bc` holds information about boundary conditions. This informationis connected to information about the triangulation, the *mesh*.Assuming $u=0$ on the boundary, we can in FEniCS generate a triangularmesh over a rectangular domain $[-1,-1]\times [-1,1]$ as follows:<jupyter_code># DO NOT RUN THIS CELL. RUN THE FULL PROGRAM (IN THE END) INSTEAD.
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 10, 10)
bc = DirichletBC(V, 0, 'on_boundary')<jupyter_output><empty_output><jupyter_text>Mathematically, the finite element method transforms our PDE toa sparse linear system. The `solve` step performs two tasks:construction of the linear system based on the given information aboutthe domain and its elements, and then solution of the linear system byeither an iterative or direct method.We are now in a position to summarize all the parts of a FEniCS programthat solves the Poisson equation by the finite element method:<jupyter_code>from fenics import *
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 10, 10)
V = FunctionSpace(mesh, 'P', 2)  # quadratic polynomials
bc = DirichletBC(V, 0, 'on_boundary')
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
u = Function(V)  # unknown FEM function to be computed
solve(a == L, u, bc)
vtkfile = File('poisson.pvd'); vtkfile << u  # store solution<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08383616748948619, 'fraction_numerical': 0.006399707441945511, 'mean_word_length': 5.218874360432063, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 26, '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': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '133036', 'license': 'MIT', 'n_tokens_mistral': 3113, 'n_tokens_neox': 2925, 'n_words': 1456, 'path': '1- overview.ipynb', 'repo_name': 'mbarzegary/finite-element-intro', 'size': '10938'}