Agglomeration example

agglomeration_poisson/CMakeLists.txt

@@ -0,0 +1,33 @@
+# Set the name of the project and target:
+SET(TARGET "agglomeration_poisson")
+
+# Declare all source files the target consists of.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  source/agglomeration_handler.cc
+  source/mapping_box.cc
+)
+
+# Add include directory
+INCLUDE_DIRECTORIES(include)
+
+
+CMAKE_MINIMUM_REQUIRED(VERSION 3.13.4)
+
+FIND_PACKAGE(deal.II 9.6.0
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
+
+# Define SOURCE_DIR for mesh file paths
+target_compile_definitions(${TARGET} PRIVATE SOURCE_DIR="${CMAKE_SOURCE_DIR}")

agglomeration_poisson/README.md

@@ -0,0 +1,135 @@
+# A Discontinuous Galerkin solver for the Poisson problem on general polytopal meshes generated through mesh agglomeration
+
+## Running the code:
+
+As in the tutorial programs, type 
+
+`cmake -DDEAL_II_DIR=/path/to/deal.II .` 
+
+on the command line to configure the program. After that you can compile with `make` and run with either `make run` or using 
+
+`./DG_advection_reaction`
+
+on the command line. 
+
+## The problem:
+This program solves the problem, for $\Omega \in \mathbb{R^2}$
+
+@f[
+\begin{cases} b \cdot \nabla u + c u = f \qquad  \text{in } \Omega \\
+\qquad \qquad u=g \qquad \text{on } \partial_{-}\Omega \end{cases}
+@f]
+
+where $g \in L^2(\partial_{-}\Omega)$ and $\partial_{-}\Omega=\{ x \in
+\partial \Omega: b(x)\cdot n(x) <0\}$ is the inflow part of the
+boundary, with $b=(b_1,b_2) \in \mathbb{R^2}$. As we know from
+classical DG theory, we need to ensure that 
+@f[
+c(x) - \frac{1}{2}\nabla \cdot b \geq \gamma_0 >0
+@f]
+for some positive $\gamma_0$ so that we have coercivity in $L^2$ at the continuous level. Discrete coercivity is achieved by using a stronger norm which takes care of jumps, see Di Pietro and Ern [2] for details.
+
+
+## The weak formulation:
+
+As trial space we choose $V_h = \{ v_h \in L^2(\Omega): v_h \in P^1(\mathbb{T_h})\} \notin H^1(\Omega)$. If we integrate by parts and sum over all cells
+
+@f[
+\sum_{T \in \mathbb{T}_h} \Bigl( (-u,\beta \cdot \nabla v_h) _T + (c
+u,v_h)_T + \bigl<(b \cdot n) u ,v_h \bigr>_{\partial T} \Bigr) =
+(f,v_h)_{\Omega}
+@f]
+
+and use the so-called DG magic formula and exploit the property $[bu]_{\mathbb{F}^i} = 0$ where $\mathbb{F}^i$ are set of internal faces we obtain the (unstable!) formulation:
+
+Find $u_h \in V_h$: 
+
+@f[
+    a_h(u_h,v_h) + b_h(u_h,v_h)=l(v_h) \qquad \forall v_h \in V_h
+@f]
+where
+@f[
+a_h(u,v_h)=\sum_{T \in \mathbb{T}_h} \Bigl( (-u,b \cdot \nabla v_h) _T + (c u,v_h)_T \Bigr)
+@f]
+
+@f[    
+b_h(u,v_h)= \sum_{F \not \in \partial_{-}\Omega} \bigl< \{ bu\}, [v_h]\bigr>_F 
+@f]
+
+@f[
+    l(v_h)= (f,v_h)_{\Omega} - \sum_{F \in \partial_{-}\Omega} \bigl< (b \cdot n) g,v_h \bigr>_F
+@f]
+
+It's well known this formulation is coercive only in $L^2$, hence the formulation is unstable as we don't "see" the derivatives. To stabilize this, we can use a jump-penalty term, i.e. our $b_h$ is replaced by:
+
+@f[
+b_h^s(u_h,v_h)=b_h(u_h,v_h)+ \sum_{F \in \mathbb{F}^i} \bigl< c_F
+[u_h],[v_h]  \bigr> 
+@f]
+
+where $c_F>0$ is a function on each edge such that $c_F \geq \theta |b \cdot n|$ for some positive $\theta$. In this program, $\theta=\frac{1}{2}$ and $c_F = \frac{1}{2} |b \cdot n|$, which corresponds to an upwind formulation. Notice that consistency is trivially achieved, as $[u]_{\mathbb{F}^i} =0$. This formulation is stable in the energy norm 
+
+@f[
+    |||\cdot ||| = \Bigl(||\cdot||_{0,\Omega}^2 + \sum_{F \in
+    \mathbb{F}}||c_F^{\frac{1}{2}}[\cdot] ||_{0,F}^2
+    \Bigr)^{\frac{1}{2}}
+@f]
+
+(well defined on $H^1(\Omega) + V_h$) and moreover we have the a-priori bound:
+
+@f[
+|||u-u_h||| \leq C h^{k+\frac{1}{2}}||u||_{k+1,\Omega} 
+@f]
+
+valid for $u \in H^{k+1}(\Omega)$.
+
+See Brezzi-Marini-Süli [3] for more details.
+
+## A-posteriori error estimator:
+
+The estimator is the one proposed by Georgoulis, Edward Hall and Charalambos Makridakis in [3]. This approach is quite different with respect to other works in the field, as the authors are trying to develop an estimator for the original hyperbolic problem, rather than taking the hyperbolic regime as the vanishing diffusivity limit.
+
+The reliability is:
+
+@f[
+|||u-u_h|||^2 \leq  C || \sqrt{b \cdot n}[u_h]||_{\Gamma^{-}}^2 + C
+\sum_{T \in \mathbb{T}_h}\Bigl( ||\beta (g-u_h^+)||_{\partial_{-}T
+\cap \partial_{-} \Omega}^2 +||f-c u_h - \Pi(f- cu_h)||_T^2 \Bigr)
+@f]
+
+where:
+
+- $\Pi$ is the (local) $L^2$ orthogonal projection onto $V_h$
+
+- $\Gamma$ is the skeleton of the mesh
+
+- $c$ is constant
+
+- $\beta = |b \cdot n|$
+
+- $u_h^+$ is the interior trace from the current cell $T$ of a the finite element function $u_h$.
+
+## Test case:
+
+The following test case has been taken from [3]. Consider:
+- $c=1$ 
+- $b=(1,1)$ 
+- $f$ to be such that the exact solution is $u(x,y)=\tanh(100(x+y-\frac{1}{2}))$
+This solution has an internal layer along the line $y=\frac{1}{2} -x$, hence we would like to see that part of the domain to be much more refined than the rest.
+
+The next image is the 3D view of the numerical solution:
+
+![Screenshot](./doc/images/warp_by_scalar_solution_layer.png)
+
+More interestingly, we see that the estimator has been able to capture the layer. Here a bulk-chasing criterion is used, with bottom fraction ´0.5´ and no coarsening. This mesh is obtained after 12 refinement cycles.
+![Screenshot](./doc/images/refined_mesh_internal_layer.png)
+
+If we look at the decrease of the energy norm of the error in the globally refined case and in the adaptively case, with respect to the DoFs, we obtain:
+
+![Screenshot](./doc/images/adaptive_vs_global_refinement.png)
+
+## References 
+* [1] Emmanuil H. Georgoulis, Edward Hall and Charalambos Makridakis (2013), Error Control for Discontinuous Galerkin Methods for First Order Hyperbolic Problems. DOI: [10.1007/978-3-319-01818-8_8
+](https://link.springer.com/chapter/10.1007%2F978-3-319-01818-8_8)
+* [2] Di Pietro, Daniele Antonio and Ern, Alexandre (2012), Mathematical Aspects of Discontinuous Galerkin Methods. ISBN: [978-3-642-22980-0](https://www.springer.com/gp/book/9783642229794)
+* [3] Franco Brezzi, Luisa Donatella Marini and Endre Süli (2004) Discontinuous Galerkin Methods for First-Order Hyperbolic Problems. DOI: [10.1142/S0218202504003866](https://doi.org/10.1142/S0218202504003866)