Galerkin FEM for elliptic PDEs
Вы не можете выбрать более 25 тем Темы должны начинаться с буквы или цифры, могут содержать дефисы(-) и должны содержать не более 35 символов.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269
  1. \chapter{Derivation}
  2. \section{Problem Statement}
  3. FEM4EllipticPDE is a Lemma module implementing a Galerkin finite element solution to a family of commonly encountered
  4. elliptic problems taking the form
  5. %Considering the uniformly elliptic Dirichlet BVP:
  6. \begin{eqnarray}
  7. \label{eq:scheme}
  8. -\nabla \cdot \left( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) \right) = g(\mathbf{r}) & \mathop{\forall}_{u \in \Omega} \\
  9. u = 0 & \mathop{\forall}_{u \in \partial \Omega}
  10. \end{eqnarray}
  11. Where $0 < \sigma_{min} \leq \sigma(\mathbf{r}) \leq \sigma_{max} \ll \infty $ and $\Omega \in \mathbb{R}^3$. When $\sigma \equiv 1$, this
  12. reduces to Poisson's equation $\left( \nabla^2 u = g \right)$. These types of problem arises in several areas of geophysics: gravitational or magnetic potentials
  13. obey Poisson's equation, while the more general form can be used to solve electrostatic (DC and SP) problems.
  14. \begin{enumerate}
  15. \item The Galerkin finite element method is defined by reposing the global problem in terms of numerous local ones through the use of
  16. appropriate test functions.
  17. \begin{itemize}
  18. \item Test functions $v$ must be from a suitable subspace.
  19. In this case, due to the Dirichlet boundary conditions,
  20. \[
  21. v \in H^a_b(\Omega)
  22. \]
  23. meet our requirements. $H^a_b$ represents the subspace of weakly differentiable functions that are zero valued at $a$, and $b$
  24. and $\textgoth{V}$ represents the actual functions. % Hilbert space.
  25. The 3D variational problem may be written
  26. % \[
  27. % -\nabla \cdot \left( \sigma(\mathbf{r}) \nabla u(\mathbf{r}) \right) = g(\mathbf{r}) & \mathop{\forall}_{u \in \Omega} \\
  28. % u = 0 & \mathop{\forall}_{u \in \partial \Omega}
  29. %
  30. %\]
  31. \[
  32. \int -v \nabla \cdot \sigma \nabla u = v g.
  33. \]
  34. This can be simplified
  35. \[
  36. \int_a^b \sigma \frac{du}{dx} \frac{dv}{dx} dx = v g
  37. \]
  38. The proof of which follows below.
  39. %\item The variational form may be found by multiplying the elliptic PDE by any arbitrary
  40. %\[
  41. % v \in H_a^b(\Omega)
  42. %\]
  43. %In 1D ($\Omega \in \mathbb{R}^1$) \autoref{eq:scheme} reduces to:
  44. %\begin{equation} \label{eq:oned}
  45. % - \DD{x} \sigma \DD{x} u = g
  46. %\end{equation}
  47. %Using the variational (weak) formulation for any arbitrary $v \in \textgoth{V}
  48. \item In any arbitrary dimension $x$, if the boundaries in that dimension of $\partial \Omega = [a,b]$ the variational problem using
  49. test functions can be constructed
  50. \begin{equation} \label{eq:oned}
  51. -\int_a^b v(x) \DD{x} \sigma(x) \DD{x} u(x) dx = \int_a^b v(x) g(x) dx ,
  52. \end{equation}
  53. The variation problem can be reduced to a simplified version, if appropriate test functions we may instead write the
  54. variational problem as
  55. \[
  56. \int_a^b \sigma \frac{du}{dx} \frac{dv}{dx} dx = \int_a^b v(x) g(x) dx
  57. \]
  58. In one dimension this resolves to
  59. \[
  60. % -\int_a^b v \frac{d}{dx} \left[ \sigma \frac{du}{dx} \right] dx =
  61. \int_a^b \sigma \frac{du}{dx} \frac{dv}{dx} dx = \int_a^b v(x) g(x) dx
  62. \]
  63. \begin{corollary}
  64. \[
  65. \frac{d}{dx} \left[ v \sigma \frac{du}{dx} \right] = v \frac{d}{dx} \left[ \sigma \frac{du}{dx} \right] +
  66. \sigma \frac{du}{dx} \frac{dv}{dx}
  67. \]
  68. \end{corollary}
  69. \begin{proof}
  70. and the integration by parts on the left hand side yields
  71. \begin{eqnarray*}
  72. \int_a^b v \frac{d}{dx} \left[ v \sigma \frac{du}{dx} \right] dx &=& \int_a^b v \frac{d}{dx}
  73. \left[ \sigma \frac{du}{dx} \right] dx +
  74. \int_a^b \sigma \frac{du}{dx} \frac{dv}{dx} dx \\
  75. \Rightarrow & & \\
  76. \left[ v \sigma \frac{du}{dx} \right]_a^b &=& \int_a^b v \frac{d}{dx} \left[
  77. \sigma \frac{du}{dx} \right] + \int_a^b \sigma \frac{dv}{dx}\frac{du}{dx} dx \\
  78. \Rightarrow & & \\\
  79. -\int_a^b v \frac{d}{dx} \left[ \sigma \frac{du}{dx} \right] dx &=& \left[ v \sigma
  80. \frac{du}{dx} \right]_a^b + \int_a^b \sigma \frac{du}{dx} \frac{dv}{dx} dx
  81. \end{eqnarray*}
  82. \[
  83. -\int_a^b v \frac{d}{dx} \left[ \sigma \frac{du}{dx} \right] dx = \int_a^b
  84. \sigma \frac{du}{dx} \frac{dv}{dx} dx
  85. \]
  86. \end{proof}
  87. Because $v\in H^1_0$ the solution vanishes at the boundaries $\left[ v \sigma \frac{du}{dx} \right]_a^b =0 $.
  88. Therefore the scheme reduces to:
  89. % Which can be rewritten as
  90. % \[
  91. % a(u,v) = <g,v>
  92. % \]
  93. The Galerkin FEM method uses triangle (hat) functions for $v$, which are renamed $\phi$ as
  94. they are specific. I will occasionally use these interchangeably.
  95. Through induction it can be shown then that the 3D problem can similarity be posed
  96. So that Equation \ref{eq:scheme} may be rewritten
  97. \[
  98. -v \nabla \cdot \sigma \nabla u = g v.
  99. \]
  100. \FloatBarrier
  101. \item \FloatBarrier Define mesh. Since $\Omega = (0,\pi)$ is an uncountable subset of $\mathbb{R}^1$ a
  102. countable subset of $\Omega$ must be defined. A mesh is defined over the interval $(0,\pi)$.
  103. This mesh is defined and shown in Figure \ref{fig:mesh}.
  104. \begin{figure}[ht]
  105. \begin{center}
  106. \input{chapters/mesh}
  107. \caption{A 1D FEM mesh is defined. $N$ nodes ($n$) are defined over the interval. Between the
  108. nodes an element is defined. The mesh is general and the spacing between nodes is defined by
  109. $N-1$ discritisation parameters $h$.}
  110. \label{fig:mesh}
  111. \end{center}
  112. \end{figure}
  113. \item Define a finite dimensional subspace based on the mesh and a particular test function.
  114. A subspace $\textgoth{V}_h \subset \textgoth{V}$ is sought such that $ \textgoth{V}_h$ is finite
  115. dimensional.
  116. Particular linearity independent test functions $\phi$ are defined on $\textgoth{V}_h$ such that
  117. $\textgoth{V}_h = \mathrm{span}(\phi_i)$. The test functions form a basis of this space.
  118. If the weak formulation of the derivative is used these
  119. test functions may be constructed such as they are only piecewise differentiable. We may define test
  120. functions as Kronecker delta $\delta$ functions taking values of unity at a particular node of the mesh.
  121. The linear spline Galerkin method of FEM interpolates these $\delta$ functions. They are zero valued at the
  122. neighbouring nodes. In 1D this interpolation yields hat functions.
  123. Therefore the test functions $\phi_i$ are defined as:
  124. \begin{equation}
  125. \phi_i(x) = \left\{ \begin{array}{ccl}
  126. 0 & \textrm{if} & x\in[\alpha, x_{i-1}] \\
  127. \frac{x - x_{i-1}}{h_i} & \textrm{if} & x\in[x_{i-1}, x_i] \\
  128. \frac{x_{i+1}-x}{h_{i+1}} & \textrm{if} & x\in[x_i, x_{i+1}] \\
  129. 0 & \textrm{if} & x\in[x_{i+1}, \beta ] \\
  130. \end{array} \right.
  131. \end{equation}
  132. Define the variational form of the solution solution
  133. \begin{equation} \label{eq:sol1}
  134. a(u, v) = <g,v>
  135. \end{equation}
  136. A is a symmetric bilinear form such that $a(u,v) = a(v,u) \forall v \in \textgoth{V}_h$.
  137. \end{itemize}
  138. In 3D things are a little different. For a given element (tetrahedra) comprised of 4 points we construct the
  139. location matrix
  140. \[ \mathbf{C} =
  141. \left[ \begin{matrix}
  142. 1 & x_1 & y_1 & z_1 \\
  143. 1 & x_2 & y_2 & z_2 \\
  144. 1 & x_3 & y_3 & z_3 \\
  145. 1 & x_4 & y_4 & z_4
  146. \end{matrix} \right]
  147. \]
  148. We can then fill the stiffness matrix
  149. \[
  150. K_{ij} = \sum_{i=0}^{N} \int_{Tet_k} = c \nabla \phi_i \cdot \phi_j
  151. \]
  152. \item The solution may be written in matrix(stiffness matrix) vector (load vector) notation.
  153. We approximate the solution (Equation \ref{eq:sol1}) using the trial and test functions $\phi$ discussed
  154. previously. Let
  155. \begin{eqnarray*}
  156. \sigma(x) &=& x + 1 \\
  157. u(x) &=& \sin(x)
  158. \end{eqnarray*}
  159. At each node $u_h \in \textgoth{V}_h : a(u_h, v) = <g,v>$ the stiffness matrix $A$
  160. \begin{eqnarray*}
  161. A_{ij} &=& a(\phi_i, \phi_j) = \int_0^\pi \sigma(x) \phi'_i \phi'_j \\
  162. &=& \int_0^\pi (x+1) \phi'_i \phi'_j
  163. \end{eqnarray*}
  164. %//Since $\phi$ is a triangle, $\phi'$ becomes a step function.
  165. We have for $i= 1,\cdots , n$
  166. \[
  167. \phi'(x) = \left\{ \begin{array}{ccl}
  168. 0 & \textrm{if} & x\in[\alpha, x_{i-1}] \\
  169. \frac{1}{h_i} & \textrm{if} & x\in[x_{i-1}, x_i] \\
  170. -\frac{1}{h_{i+1}} & \textrm{if} & x\in[x_i, x_{i+1}] \\
  171. 0 & \textrm{if} & x\in[x_{i+1}, \beta ] \\
  172. \end{array} \right.
  173. \]
  174. The product $\phi_i'(x) \phi_j(x)$ vanishes when $|i - j| > 1$. So that in a 1D case, A becomes tri-diagonal.
  175. Following the general equation:
  176. \begin{eqnarray*}
  177. A_{ij} &=& \int_{i-1}^{i+1} (x+1) \phi'_i(x) \phi'_j(x) \\
  178. \end{eqnarray*}
  179. The load vector $g$ is given:
  180. \[
  181. [g]_i = <g, \phi_i> = \int_0^\pi g \phi_i
  182. \]
  183. Using $ g = [g]$, the linear system becomes:
  184. \[
  185. Au = g
  186. \]
  187. \item For this case the entries of the stiffness matrix and load vector are
  188. For the three non-zero cases of $A$:
  189. \begin{eqnarray*}
  190. A_{ii} &=& \int_{i-1}^{i+1} (x+1) \phi'_i(x) \phi'_i(x) \\
  191. A_{ii} &=& \int_{{i-1}}^{i} (x+1) \phi_i'^2 dx + \int_{i}^{x_{i+1}} (x+1)\phi_i'^2 \\
  192. A_{ii} &=& \frac{1}{h_i^2}\int_{{i-1}}^{i} (x+1) dx + \frac{1}{h_{i+1}^2} \int_{i}^{{i+1}} (x+1) dx
  193. \end{eqnarray*}
  194. \begin{eqnarray*}
  195. A_{i,i+1} &=& \int_{i}^{i+1} (x+1) \phi'_i(x) \phi'_{i+1}(x) \\
  196. A_{i,i+1} &=& \int_{i}^{i+1} (x+1) \phi_i' \phi_{i+1}' dx \\ %+ \int_{i}^{x_{i+1}} (x+1)\phi_i' \phi_{i+1}' \\
  197. A_{i,i+1} &=& \frac{1}{h_i + h_{i+1}}\int_{i}^{i+1} (x+1) dx %+ \frac{1}{h_{i+1}^2} \int_{i}^{{i+1}} (x+1) dx
  198. \end{eqnarray*}
  199. \begin{eqnarray*}
  200. A_{i,i-1} &=& \int_{i-1}^{i} (x+1) \phi'_i(x) \phi'_{i-1}(x) \\
  201. A_{i,i-1} &=& \int_{i-1}^{i} (x+1) \phi_i' \phi_{i-1}' dx \\%+ \int_{i}^{x_{i+1}} (x+1)\phi_i' \phi_{i-1}' \\
  202. A_{i,i-1} &=& \frac{1}{h_i + h_{i-1}}\int_{{i-1}}^{i} (x+1) dx % dx + \frac{1}{h_{i+1}^2} \int_{i}^{{i+1}} (x+1) dx
  203. \end{eqnarray*}
  204. These integrals are evaluated numerically in my program \textit{gfemddn} using Simpson's Rule.
  205. Solving for the forcing function $g$ in \autoref{eq:scheme}
  206. \begin{eqnarray*}
  207. -\nabla \cdot \sigma(x) \nabla u(x) &=& g(x) \\
  208. -\frac{\partial}{\partial x} \cdot (x+1) \frac{\partial }{\partial x} \sin(x) &=& g(x) \\
  209. -\frac{\partial}{\partial x} \cdot (x+1) \cos(x) &=& g(x) \\
  210. \sin(x) + x\sin(x) - \cos(x) &=& g(x) \\
  211. \sin(x) (x+1) - \cos(x) &=& g(x) \\
  212. \end{eqnarray*}
  213. The load vector $g$ is then:
  214. \[
  215. g_i = <g, \phi_i> = \int_0^\pi \left( \sin(x) (x+1) - \cos(x) \right) \phi_i
  216. \]
  217. In \textit{gfemddn} this integral is again evaluated using Simpson's Rule.
  218. \end{enumerate}