You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
# - How to solve a simple with a Discontinuous Galerkin (DG) discretization
9
+
# - How to build discontinuous FE spaces
10
+
# - How to integrate quantities on the mesh skeleton
11
+
# - How to compute jumps and averages of quantities on the mesh skeleton
12
+
#
13
+
#
14
+
# ## Problem statement
15
+
#
16
+
# The goal of this tutorial is to solve a simple PDE using a Discontinuous Galerkin (DG) formulation. For simplicity, we take the Poisson equation on the unit multi dimensional cube $\Omega \doteq (0,1)^d$, with $d=2$ and $d=3$, as a model problem:
17
+
#
18
+
#
19
+
# ```math
20
+
# \left\lbrace
21
+
# \begin{aligned}
22
+
# -\Delta u = f \ \text{in} \ \Omega\\
23
+
# u = g \ \text{on}\ \partial\Omega,\\
24
+
# \end{aligned}
25
+
# \right.
26
+
# ```
27
+
# where $f$ is the source term and $g$ is the Dirichlet boundary value.
28
+
#
29
+
# We are going to solve two version of this problem. In a first stage, we take $d=3$ and consider a manufactured solution, namely $u(x) = 3 x_1 + x_2 + 2 x_3$, that belongs to the FE interpolation that we will build below. In this case, we expect to compute a numerical solution with an approximation error close to the machine precision. On the other hand, we will perform a convergence test for the 2D case ($d=2$). To this, end we will consider a manufactured solution that cannot by represented exactly by the interpolation, namely $u(x)=x_2 \sin(2 \pi\ x_1)$. Our goal is to confirm that the convergence order of the discretization error is the optimal one.
30
+
#
31
+
# ## Numerical Scheme
32
+
#
33
+
# In contrast to previous tutorials, we consider a DG formulation to approximate the problem. For the sake of simplicity, we take the well know (symmetric) interior penalty method. For this formulation, the approximation space is made of discontinuous piece-wise polynomials, namely
34
+
#
35
+
# ```math
36
+
# V \doteq \{ v\in L^2(\Omega):\ v|_{T}\in Q_p(T) \text{ for all } T\in\mathcal{T} \},
37
+
# ```
38
+
# where $\mathcal{T}$ is the set of all cells $T$ of the FE mesh, and $Q_p(T)$ is a polynomial space of degree $p$ defined on a generic cell $T$. For simplicity, we consider Cartesian meshes in this tutorial. In this case, the space $Q_p(T)$ is made of multi-variate polynomials up to degree $p$ in each spatial coordinate.
39
+
#
40
+
#
41
+
#
42
+
#
43
+
# In order to write the weak form of the problem, we need to introduce the set of interior and boundary facets associated with the FE mesh, denoted here as $\mathcal{F}_\Gamma$ and $\mathcal{F}_{\partial\Omega}$ respectively. In addition, for a given function $v\in V$ restricted to the interior facets $\mathcal{F}_\Gamma$, we need to define the well known jump and mean value operators:
# with $v^+$, and $v^-$ being the restrictions of $v\in V$ to the cells $T^+$, $T^-$ that share a generic interior facet in $\mathcal{F}_\Gamma$, and $n^+$, and $n^-$ are the facet outward unit normals from either the perspective of $T^+$ and $T^-$ respectively.
48
+
#
49
+
# With this notation, the weak form associated with the interior penalty formulation reads: find $u\in V$ such that $a(v,u) = b(v)$ for all $v\in V$. The bilinear $a(\cdot,\cdot)$ and linear form $b(\cdot)$ have contributions associated with the bulk of $\Omega$, and the boundary and interior facets $\mathcal{F}_{\partial\Omega}$, $\mathcal{F}_\Gamma$, namely
# for the interior facets. In previous expressions, $|F|$ denotes the diameter of the face $F$ (in our Cartesian grid, this is equivalent to the characteristic mesh size $h$), and $\gamma$ is a stabilization parameter that should be chosen large enough such that the bilinear form $a(\cdot,\cdot)$ is stable and continuous. Here, we take $\gamma = p\ (p+1)$.
72
+
#
73
+
# ## 3D manufactured solution
1
74
2
75
using Gridap
3
76
import Gridap: ∇
4
77
5
-
u(x) = x[1] + x[2]
6
-
∇u(x) =VectorValue(1.0,1.0)
78
+
u(x) =3*x[1] + x[2] +2*x[3]
79
+
∇u(x) =VectorValue(3.0,1.0,2.0)
7
80
∇(::typeof(u)) = ∇u
8
81
f(x) =0.0
9
82
g(x) =u(x)
10
83
11
84
L =1.0
12
-
limits = (0.0, L, 0.0, L)
13
-
ncellx=4
14
-
model =CartesianDiscreteModel(domain=limits, partition=(ncellx,ncellx))
85
+
limits = (0.0, L, 0.0, L, 0.0, L)
86
+
n=4
87
+
model =CartesianDiscreteModel(domain=limits, partition=(n,n,n))
0 commit comments