The aim of this test is:
We solve Laplace’s Equation in 2D on a $1 \times 1$ square domain. The boundary conditions will be specified below.
Laplace’s equation is
$$ \begin{equation} \ \mathop{\mathrm{div}} (a \mathop{\mathrm{grad}} u) = 0 \end{equation} $$The weak form is derived as usual by multiplying with a test function $v$ and integrating over the domain $\Omega$:
$$ \begin{equation} \ \int_{\Omega} v \mathop{\mathrm{div}} (a \mathop{\mathrm{grad}} u) \, \mathrm{d}\Omega = 0 \,, \end{equation} $$which can be transformed further to
$$ \begin{equation} \int_{\Omega} a \mathop{\mathrm{grad}} v \cdot \mathop{\mathrm{grad}} u \, \mathrm{d}\Omega = \int_{\Omega} \mathop{\mathrm{div}} (v a \mathop{\mathrm{grad}} u) \, \mathrm{d}\Omega = \int_{\Gamma_{\mathrm{N}}} v a \mathop{\mathrm{grad}} u \cdot n \, \mathrm{d}\Gamma \,, \end{equation} $$where in the second equality Gauss’s theorem has been applied. As usual, the domain boundary $\partial\Omega = \Gamma\_{\mathrm{D}} \cup \Gamma\_{\mathrm{N}}$ is subdivided into the Dirichlet and the Neumann boundary and $v$ vanishes on $\Gamma\_{\mathrm{D}}$. The r.h.s. of the above equation is the total flux associated with $u$ flowing into the domain $\Omega$ through $\Gamma\_{\mathrm{N}}$: $a \mathop{\mathrm{grad}} u$ is the flux density and $\n$ is the inwards directed surface normal.
The weak form just derived is implemented (after FEM discretization) in the groundwater flow process in OpenGeoSys. Note that for the application of Neumann boundary conditions, it is necessary to know whether the flux has to be applied with a positive or a negative sign!
The coefficient $a$ of Laplace’s equation is taken to be unity. By differentiation it can be easily checked that
$$ \begin{equation} u(x, y) = \sin(bx) \sinh(by) \end{equation} $$solves Laplace’s equation inside $\Omega$ for any $b$. In this example we set $b = \tfrac 23 \pi$.
As boundary conditions we apply Dirichlet BCs at the top, left and bottom of the domain with values from $u(x,y)\_{\Gamma_{\mathrm{D}}}$. On the right boundary of the domain a Neumann BC is applied. There $n = (1, 0)$, which implies that $a \mathop{\mathrm{grad}} u \cdot n = a \, \partial u / \partial x$.
The numerical result obtained from OpenGeoSys is:
The absolute difference between the analytical and numerical solutions is smaller than $4 \cdot 10^{4}$:
This article was written by Christoph Lehmann. If you are missing something or you find an error please let us know.
Generated with Hugo 0.122.0
in CI job 437561

Last revision: April 3, 2024
Commit: [PL/CT] Pass shape matrix cache to local assemblers of ComponentTransport 845cbb1
 Edit this page on