* Constraints typically come from several sources, for example:
* - If you have Dirichlet-type boundary conditions, $u|_{\partial\Omega}=g$,
* one usually enforces
- * them by requiring that that degrees of freedom on the boundary have
+ * them by requiring that degrees of freedom on the boundary have
* particular values, for example $x_{12}=42$ if the boundary condition
* $g(\mathbf x)$ requires that the finite element solution $u(\mathbf x)$
* at the location of degree
* below. The resulting linear systems are equivalent to those one gets after
* calling the ConstraintMatrix::condense() functions.
*
+ * @note Both ways of applying constraints set the value of the matrix
+ * diagonals to constrained entries to a <i>positive</i> entry of the same
+ * magnitude as the other entries in the matrix. As a consequence, you need to
+ * set up your problem such that the weak form describing the main matrix
+ * contribution is not <i>negative definite</i>. Otherwise, iterative solvers
+ * such as CG will break down or be considerably slower as GMRES.
+ *
* @note While these two ways are <i>equivalent</i>, i.e., the solution of
* linear systems computed via either approach is the same, the linear
* systems themselves do not necessarily have the same matrix and right
* as a result of the way in which we compute them; they are, however, always
* chosen in such a way that the solution to the linear system is the same.
*
- *
* <h4>Condensing matrices and sparsity patterns</h4>
*
* As mentioned above, the first way of using constraints is to build linear
<h3>Specific improvements</h3>
<ol>
+
+ <li>
+ Changed: find_active_cell_around_point now throws the exception ExcPointNotFound
+ if the point is outside the Triangulation. This exception can be caught.
+ <br>
+ (Timo Heister, 2013/09/18)
+ </li>
+
+ <li>
+ Changed: we now call MPI_Init_thread instead of MPI_Init.
+ <br>
+ (Timo Heister, 2013/09/17)
+ </li>
+
<li>
Enhancement: It is now possible to use the build directory directly without
the need to install first. For this, a second copy of all necessary project
submitted, 2013.
</li>
+ <li> G. Alam
+ <br>
+ <strong>Fast iterative solution of large scale statistical
+ inverse problems
+ </strong>
+ <br>
+ Master Thesis, Stockholm University, 2013.
+ </li>
+
<li> J. C. Araujo-Cabarcas
<br>
<strong>Numerical methods for glacial isostatic adjustment models
MSc thesis, Uppsala University, Uppsala, Sweden, 2013.
</li>
- <li> G. Alam
+ <li> P. Argul, G. Ruocci, F. Freddi, K. Benzarti
<br>
- <strong>Fast iterative solution of large scale statistical
- inverse problems
+ <strong>An improved damage modelling to deal with the variability of
+ fracture mechanisms in FRP reinforced concrete structures
</strong>
<br>
- Master Thesis, Stockholm University, 2013.
+ International Journal of Adhesion and Adhesives 45 , 7-20, 2013.
</li>
<li> D. Arndt
Msc thesis, University of Göttingen, Germany, 2013.
</li>
+ <li> A. Barker, T. Rees, M. Stoll
+ <br>
+ <strong>A fast solver for an H<sub>1</sub> regularized
+ PDE-constrained optimization problem</strong>
+ <br>
+ submitted, 2013.
+ </li>
+
+ <li> S. Bartels, R. H. Nochetto, A. J. Salgado
+ <br>
+ <strong>A total variation diminishing interpolation operator
+ and applications</strong>
+ <br>
+ submitted, 2013.
+ </li>
+
<li> P. Benner, M. Stoll
<br>
<strong>Optimal control for Allen-Cahn equations enhanced by model predictive control</strong>
Computers & Mathematics with Applications (in press), 2013.
</li>
+ <li> F. Freddi, G. Royer-Carfagni
+ <br>
+ <strong>Plastic flow as an energy minimization problem. Numerical
+ Experiments
+ </strong>
+ <br>
+ Journal of Elasticity, in press, 2013.
+ </li>
+
<li> A. Geringer, B. Lenhof, S. Diebels
<br>
<strong>Macroscopic modeling of foams: an order-parameter
pp. 46-72, 2013.
</li>
+ <li> S. Kim
+ <br>
+ <strong>Analysis of the convected Helmholtz equation with a
+ uniform mean flow in a waveguide with complete radiation
+ boundary conditions
+ </strong>
+ <br>
+ Journal of Mathematical Analysis and Applications,
+ in press,
+ <!-- http://www.sciencedirect.com/science/article/pii/S0022247X1300749X -->
+ 2013.
+ </li>
+
<li> N. Klein
<br>
<strong>Consistent FE-Analysis of elliptic Variational Inequalities
MSc Thesis, University of Cape Town, 2013.
</li>
+ <li> T. Reis, W. Wollner
+ <br>
+ <strong>Iterative solution of operator Lyapunov equations
+ arising in heat transfer
+ </strong>
+ <br>
+ 2013 European Control Conference, Zurich Switzerland,
+ pp. 41-46, 2013.
+ </li>
+
<li> T. Richter, T. Wick
<br>
<strong>Optimal Control and Parameter Estimation for Stationary Fluid-Structure Interaction Problems
must solve in an implicit system. This is because, unlike continuous finite
elements, in typical discontinuous elements there is one degree of freedom at
each vertex <i>for each of the adjacent elements</i>, rather than just one,
-and similarly for edges and faces. As another example,
-for the FE_DGP_Monomial basis, each
+and similarly for edges and faces. As an example of how fast the number of
+unknowns grows,
+consider the <code>FE_DGP_Monomial</code> basis: each
scalar solution component is represented by polynomials of degree $p$
-which yields $(1/dim!)*\prod_{i=1}^{dim}(p+i)$ degrees of freedom per
-element. Typically, all degrees of freedom in an element are coupled
+with $(1/dim!)*\prod_{i=1}^{dim}(p+i)$ degrees of freedom per
+element. Typically, all degrees of freedom in an element are coupled
to all of the degrees of freedom in the adjacent elements. The resulting
discrete equations yield very large linear systems very quickly, especially
-for systems of equations in dim=2 or dim=3.
+for systems of equations in 2 or 3 dimensions.
<h4> Reducing the size of the linear system </h4>
-To alleviate the computational cost of solving such large linear systems,
+To alleviate the computational cost of solving such large linear systems,
the hybridizable discontinuous Galerkin (HDG) methodology was introduced
-by Cockburn and co-workers
-(see the references in the recent HDG overview article by
+by Cockburn and co-workers
+(see the references in the recent HDG overview article by
N.C. Nguyen and J. Peraire:
- <i>Hybridizable discontinuous Galerkin methods for partial differential
+ <i>Hybridizable discontinuous Galerkin methods for partial differential
equations in continuum mechanics</i>, Journal of Computational Physics,
2012, 231:18, 5955-5988.
<a href="http://dx.doi.org/10.1016/j.jcp.2012.02.033">[DOI]</a>).
mappings. The partial differential equations are first written as a first
order system, and each field is then discretized via a DG method. At this
point the single-valued "trace" values on the skeleton of the
-mesh, i.e. element faces, are taken to be independent unknown quantities.
+mesh, i.e. element faces, are taken to be independent unknown quantities.
The Dirichlet-to-Neumann map concept then permits the following solution procedure:
<ol>
- <li> Use local element interior data to enforce a Neumann condition on the
-skeleton of the triangulation. The global problem is then to solve for the
+ <li> Use local element interior data to enforce a Neumann condition on the
+skeleton of the triangulation. The global problem is then to solve for the
trace values, which are the only globally coupled unknowns.
<li> Use the known skeleton values as Dirichlet data for solving local
element-level solutions. This is known as the
-'local solver', and is an <i>embarrassingly parallel</i> element-by-element
+'local solver', and is an <i>embarrassingly parallel</i> element-by-element
solution process.
</ol>
\begin{pmatrix} A & B \\ C & D \end{pmatrix}
\begin{pmatrix} U \\ \Lambda \end{pmatrix}
=
-\begin{pmatrix} F \\ G \end{pmatrix}
+\begin{pmatrix} F \\ G \end{pmatrix}.
@f}
Our aim is now to eliminate the <i>U</i> block with a Schur complement
approach similar to step-20, which results in the following two steps:
@f{eqnarray*}
-(D - C A^{-1} B) \Lambda &=& G - C A^{-1} F \\
-A U &=& F - B \Lambda
+(D - C A^{-1} B) \Lambda &=& G - C A^{-1} F, \\
+A U &=& F - B \Lambda.
@f}
The point is that the presence of $A^{-1}$ is not a problem because $A$ is a
block diagonal matrix where each block corresponds to one cell and is
matrix <i>A</i> element by element (the local solution of the Dirichlet
problem) and subtract $CA^{-1}B$ from $D$. The steps in the Dirichlet-to-Neumann map concept hence correspond to
<ol>
- <li> constructing the Schur complement matrix $D-C A^{-1} B$ and right hand side $G - C A^{-1} F$,
+ <li> constructing the Schur complement matrix $D-C A^{-1} B$ and right hand side $G - C A^{-1} F$ <i>locally on each cell</i>
+ and inserting the contribution into the global trace matrix in the usual way,
<li> solving the Schur complement system for $\Lambda$, and
- <li> solving the equation for <i>U</i> using the second equation which uses $\Lambda$.
+ <li> solving for <i>U</i> using the second equation, given $\Lambda$.
</ol>
<h4> Solution quality and rates of convergence</h4>
Another criticism of traditional DG methods is that the approximate fluxes
converge suboptimally. The local HDG solutions can be shown to converge
-as $\mathcal{O}(h^{p+1})$, i.e., at optimal order. Additionally, a
-super-convergence property can be used to post-process a new approximate
+as $\mathcal{O}(h^{p+1})$, i.e., at optimal order. Additionally, a
+super-convergence property can be used to post-process a new approximate
solution that converges at the rate $\mathcal{O}(h^{p+2})$.
<h3> HDG applied to the convection-diffusion problem </h3>
-The HDG formulation used for this example is taken from
+The HDG formulation used for this example is taken from
<br>
<b>
N.C. Nguyen, J. Peraire, B. Cockburn:
- <i>An implicit high-order hybridizable discontinuous Galerkin method
- for linear convection–diffusion equations</i>,
- Journal of Computational Physics, 2009, 228:9, 3232-3254.
+ <i>An implicit high-order hybridizable discontinuous Galerkin method
+ for linear convection–diffusion equations</i>,
+ Journal of Computational Physics, 2009, 228:9, 3232-3254.
<a href="http://dx.doi.org/10.1016/j.jcp.2009.01.030">[DOI]</a>
</b>
-We consider the convection-diffusion equation:
+We consider the convection-diffusion equation over the domain $\Omega$
+with Dirichlet boundary $\partial \Omega_D$ and Neumann boundary
+$\partial \Omega_N$:
@f{eqnarray*}
\nabla \cdot (\mathbf{c} u) - \nabla \cdot (\kappa \nabla u) &=& f,
\quad \text{ in } \Omega, \\
u &=& g_D, \quad \text{ on } \partial \Omega_D, \\
- (\mathbf{c} u - \kappa \nabla u)\cdot \mathbf{n} &=& g_N,
+ (\mathbf{c} u - \kappa \nabla u)\cdot \mathbf{n} &=& g_N,
\quad \text{ on } \partial \Omega_N.
@f}
\mathbf{q} + \kappa \nabla u &=& 0, \quad \text{ in } \Omega, \\
\nabla \cdot (\mathbf{c} u + \mathbf{q}) &=& f, \quad \text{ in } \Omega, \\
u &=& g_D, \quad \text{ on } \partial \Omega_D, \\
- (\mathbf{q} + \mathbf{c}u)\cdot\mathbf{n} &=& g_N,
+ (\mathbf{q} + \mathbf{c}u)\cdot\mathbf{n} &=& g_N,
\quad \text{ on } \partial \Omega_N.
@f}
-We multiply these equations by the weight functions $\mathbf{v}, w$
+We multiply these equations by the weight functions $\mathbf{v}, w$
and integrate by parts over every element $K$ to obtain:
@f{eqnarray*}
(\mathbf{v}, \kappa^{-1} \mathbf{q})_K - (\nabla\cdot\mathbf{v}, u)_K
+ \left<\mathbf{v}\cdot\mathbf{n}, \hat{u}\right>_{\partial K} &=& 0, \\
- - (\nabla w, \mathbf{c} u + \mathbf{q})_K
- + \left<(w, \hat{\mathbf{c} u}+\hat{\mathbf{q}})\cdot\mathbf{n}\right>_{\partial K}
+ - (\nabla w, \mathbf{c} u + \mathbf{q})_K
+ + \left<(w, \widehat{\mathbf{c} u}+\hat{\mathbf{q}})\cdot\mathbf{n}\right>_{\partial K}
&=& (w,f)_K.
@f}
The terms decorated with a hat denote the numerical traces (also commonly referred
-to as numerical fluxes). They are approximations
-to the interior values on the boundary of the element. To ensure conservation,
+to as numerical fluxes). They are approximations
+to the interior values on the boundary of the element. To ensure conservation,
these terms must be single-valued on any given element edge $\partial K$ even
though, with discontinuous shape functions, there may of course be multiple
values coming from the cells adjacent to an interface.
We eliminate the numerical trace $\hat{\mathbf{q}}$ by using traces of the form:
@f{eqnarray*}
- \hat{\mathbf{c} u}+\hat{\mathbf{q}} = \mathbf{c}\hat{u} + \mathbf{q}
+ \widehat{\mathbf{c} u}+\hat{\mathbf{q}} = \mathbf{c}\hat{u} + \mathbf{q}
+ \tau(u - \hat{u})\mathbf{n} \quad \text{ on } \partial K.
@f}
+
+The variable $\hat {u}$ is introduced as an additional independent variable
+and is the one for which we finally set up a globally coupled linear
+system. As mentioned above, it is defined on the element faces and
+discontinuous from one face to another.
+Values for $u$ and $\mathbf{q}$ appearing in the numerical trace function
+are taken to be the cell's interior solution restricted
+to the boundary $\partial K$.
+
The local stabilization parameter $\tau$ has effects on stability and accuracy
of HDG solutions; see the literature for a further discussion. A stabilization
parameter of unity is reported to be the choice which gives best results. A
stabilization parameter $\tau$ that tends to infinity prohibits jumps in the
solution over the element boundaries, making the HDG solution approach the
-approximation of continuous finite elements. In the program below, we choose
+approximation with continuous finite elements. In the program below, we choose
the stabilization parameter as
@f{eqnarray*}
\tau = \frac{\kappa}{\ell} + |\mathbf{c} \cdot \mathbf{n}|
@f}
where we set the diffusion length scale to $\ell = \frac{1}{5}$.
-The trace/skeleton variables in HDG methods are single-valued on element faces. As such,
-they must strongly represent the Dirichlet data on $\partial\Omega_D$. We introduce
-a new variable $\lambda$ such that
+The trace/skeleton variables in HDG methods are single-valued on element
+faces. As such, they must strongly represent the Dirichlet data on
+$\partial\Omega_D$. This means that
@f{equation*}
- \hat{u} =
- \begin{cases}
- g_D & \text{ on } \partial \Omega_D, \\
- \lambda & \text{otherwise}.
-\end{cases}
-@f}
+ \hat{u}|_{\partial \Omega_D} = g_D,
+@f}
+where the equal sign actually means an $L_2$ projection of the boundary
+function $g$ onto the space of the face variables (e.g. linear functions on
+the faces). This constraint is then applied to the skeleton variable $\hat{u}$
+using inhomogeneous constraints by the method
+VectorTools::project_boundary_values.
-Eliminating $\hat{u}$ from the weak form in favor of $\lambda$, summing the elemental
+Summing the elemental
contributions across all elements in the triangulation, enforcing the normal
component of the numerical flux, and integrating by parts
-on the equation weighted by $w$, we arrive at the final form of the problem:
-Find $(\mathbf{q}_h, u_h, \lambda_h) \in
+on the equation weighted by $w$, we arrive at the final form of the problem:
+Find $(\mathbf{q}_h, u_h, \hat{u}_h) \in
\mathcal{V}_h^p \times \mathcal{W}_h^p \times \mathcal{M}_h^p$ such that
@f{eqnarray*}
(\mathbf{v}, \kappa^{-1} \mathbf{q}_h)_{\mathcal{T}}
- ( \nabla\cdot\mathbf{v}, u_h)_{\mathcal{T}}
- + \left<\mathbf{v}\cdot\mathbf{n}, \lambda_h\right>_{\partial\mathcal{T}}
- &=&
- - \left<\mathbf{v}\cdot\mathbf{n}, g_D\right>_{\partial\Omega_D},
+ + \left<\mathbf{v}\cdot\mathbf{n}, \hat{u}_h\right>_{\partial\mathcal{T}}
+ &=& 0,
\quad \forall \mathbf{v} \in \mathcal{V}_h^p,
\\
(w, \mathbf{c}\nabla u_h + \nabla \cdot \mathbf{q}_h)_{\mathcal{T}}
- + \left<w, \tau (u_h - \lambda_h)\right>_{\partial \mathcal{T}}
+ + \left<w, \tau (u_h - \hat{u}_h)\right>_{\partial \mathcal{T}}
&=&
(w, f)_{\mathcal{T}},
\quad \forall w \in \mathcal{W}_h^p,
\\
- \left< \mu, \mathbf{c} \lambda_h\cdot \mathbf{n}
+ \left< \mu, \hat{u}_h\mathbf{c} \cdot \mathbf{n}
+ \mathbf{q}_h\cdot \mathbf{n}
- + \tau (u_h - \lambda_h)\right>_{\partial \mathcal{T}}
+ + \tau (u_h - \hat{u}_h)\right>_{\partial \mathcal{T}}
&=&
\left<\mu, g_N\right>_{\partial\Omega_N},
\quad \forall \mu \in \mathcal{M}_h^p.
@f}
-The unknowns $(\mathbf{q}_h, u_h)$ are referred to as local variables; they are
-represented as standard DG variables. The unknown $\lambda_h$ is the skeleton
+The unknowns $(\mathbf{q}_h, u_h)$ are referred to as local variables; they are
+represented as standard DG variables. The unknown $\hat{u}_h$ is the skeleton
variable which has support on the codimension-1 surfaces (faces) of the mesh.
+
+We use the notation $(\cdot, \cdot)_{\mathcal{T}} = \sum_K (\cdot, \cdot)_K$
+to denote the sum of integrals over all cells and $\left<\cdot,
+\cdot\right>_{\partial \mathcal{T}} = \sum_K \left<\cdot,
+\cdot\right>_{\partial K}$ to denote integration over all faces of all cells,
+i.e., interior faces are visited twice. When combining the contribution from
+both elements sharing a face, the above equation yields terms familiar for DG
+with jumps of the solution over the cell boundaries.
+
+In the equation above, the space $\mathcal {W}_h^{p}$ for the scalar variable
+<i>u<sub>h</sub></i> is defined as the space of functions that are tensor
+product polynomials of degree <i>p</i> on each cell and discontinuous over the
+element boundaries $\mathcal Q_{-p}$, i.e., the space described by
+<code>FE_DGQ<dim>(p)</code>. The space for the gradient or flux variable
+<b>q</b><i><sub>h</sub></i> is a vector element space where each component is
+a locally polynomial and discontinuous $\mathcal Q_{-p}$. In the code below,
+we collect these two local parts together in one FESystem where the first @p
+dim components denote the gradient part and the last scalar component
+corresponds to the scalar variable. For the skeleton component $\hat{u}_h$, we
+define a space that consists of discontinuous tensor product polynomials that
+live on the element faces, which in deal.II is implemented by the class
+FE_FaceQ. This space is otherwise similar to FE_DGQ, i.e., the solution
+function is not continuous between two neighboring faces, see also the results
+section below for an illustration.
+
In the weak form given above, we can note the following coupling patterns:
<ol>
<li> The matrix $A$ consists of local-local coupling terms. These arise when the
- local weighting functions $(\mathbf{v}, w)$ multiply the local solution terms
+ local weighting functions $(\mathbf{v}, w)$ multiply the local solution terms
$(\mathbf{q}_h, u_h)$.
- <li> The matrix $B$ represents the local-face coupling. These are the terms
- with weighting functions $(\mathbf{v}, w)$ multiplying the skeleton variable
- $\lambda_h$.
- <li> The matrix $C$ represents the face-local coupling, which involves the
+ <li> The matrix $B$ represents the local-face coupling. These are the terms
+ with weighting functions $(\mathbf{v}, w)$ multiplying the skeleton variable
+ $\hat{u}_h$.
+ <li> The matrix $C$ represents the face-local coupling, which involves the
weighting function $\mu$ multiplying the local solutions $(\mathbf{q}_h, u_h)$.
- <li> The matrix $D$ is the face-face coupling;
- terms involve both $\mu$ and $\lambda_h$.
+ <li> The matrix $D$ is the face-face coupling;
+ terms involve both $\mu$ and $\hat{u}_h$.
</ol>
<h4> Post-processing and super-convergence </h4>
+One special feature of the HDG methods is that they typically allow for
+constructing an enriched solution that gains accuracy. This post-processing
+takes the HDG solution in an element-by-element fashion and combines it such
+that one can get $\mathcal O(h^{p+2})$ order of accuracy when using
+polynomials of degree <i>p</i>. For this to happen, there are two necessary
+ingredients:
+<ol>
+ <li> The computed solution gradient $\mathbf{q}_h$ converges at optimal rate,
+ i.e., $\mathcal{O}(h^{p+1})$.
+ <li> The average of the scalar part of the solution, <i>u<sub>h</sub></i>,
+ on each cell <i>K</i> super-converges at rate $\mathcal{O}(h^{p+2})$.
+</ol>
+
+We now introduce a new variable $u_h^* \in \mathcal{V}_h^{p+1}$, which we find
+by the expression $|\kappa \nabla u_h^* + \mathbf{q}_h|^2$ over the cell
+<i>K</i> under the constraint $\left(1, u_h^*\right)_K &=& \left(1,
+u_h\right)_K$. This translates to the following system of equations:
+@f{eqnarray*}
+\left(1, u_h^*\right)_K &=& \left(1, u_h\right)_K\\
+\left(\nabla w_h^*, \kappa \nabla u_h^*\right)_K &=&
+-\left(\nabla w_h^*, \mathbf{q}_h\right)_K
+\quad \text{for all } w_h^* \in \mathcal Q^{p+1}.
+@f}
+
+Since we test by the whole set of basis functions in the space of tensor
+product polynomials of degree <i>p</i>+1 in the second set of equations, this
+is an overdetermined system with one more equation than unknowns. We fix this
+in the code below by omitting one of these equations (since the rows in the
+Laplacian are linearly dependent when representing a constant function). As we
+will see below, this form of the post-processing gives the desired
+super-convergence result with rate $\mathcal {O}(h^{p+2})$. It should be
+noted that there is some freedom in constructing $u_h^*$ and this minimization
+approach to extract the information from the gradient is not the only one. In
+particular, the post-processed solution defined here does not satisfy the
+convection-diffusion equation in any sense. As an alternative, the paper by
+Nguyen, Peraire and Cockburn cited above suggests another somewhat more
+involved formula for convection-diffusion that can also post-process the flux
+variable into an $H(\Omega,\mathrm{div})$-conforming variant and better
+represents the local convection-diffusion operator when the diffusion is
+small. We leave the implementation of a more sophisticated post-processing as
+a possible extension to the interested reader.
+
+Note that for vector-valued problems, the post-processing works similarly. One
+simply sets the constraint for the mean value of each vector component
+separately and uses the gradient as the main source of information.
<h3> Problem specific data </h3>
-Example specific data:
-$\Omega := [-1,1]^d$
+For this tutorial program, we consider almost the same test case as in
+step-7. The computational domain is $\Omega := [-1,1]^d$ and the exact
+solution corresponds to the one in step-7, except for a scaling. We use the
+following source centers <i>x<sub>i</sub></i> for the exponentials
+<ul>
+ <li> 1D: $\{x_i\}^1 = \{ -\frac{1}{3}, 0, \frac{1}{3} \}$,
+ <li> 2D: $\{\mathbf{x}_i\}^2 = \{ (-\frac{1}{2},\frac{1}{2}),
+ (-\frac{1}{2},-\frac{1}{2}),
+ (\frac{1}{2},-\frac{1}{2})
+ \}$,
+ <li> 3D: $\{\mathbf{x}_i\}^3 = \{ (-\frac{1}{2},\frac{1}{2}, \frac{1}{4}),
+ (-\frac{3}{5},-\frac{1}{2}, -\frac{1}{8}),
+ (\frac{1}{2},-\frac{1}{2}, \frac{1}{2})
+ \}$.
+</ul>
+
+With the exact solution given, we then choose the forcing on the right hand
+side and the Neumann boundary condition such that we obtain this solution
+(manufactured solution technique). In this example, we choose the diffusion
+equal to one and the convection as
\f[
\mathbf{c} = \begin{cases}
1, & \textrm{dim}=1 \\
(y, -x, 1), & \textrm{dim}=3
\end{cases}
\f]
-
-the centers $x_i$ of the exponentials are
- $\{x_i\}_1^3 = \{ -\frac{1}{3}, 0, \frac{1}{3} \}$
- $\{\mathbf{x}_i\}_1^3 = \{ (-\frac{1}{2},\frac{1}{2})
- (-\frac{1}{2},-\frac{1}{2})
- (\frac{1}{2},-\frac{1}{2})
- \}$
- $\{\mathbf{x}_i\}_1^3 = \{ (-\frac{1}{2},\frac{1}{2}, \frac{1}{4})
- (-\frac{3}{5},-\frac{1}{2}, -\frac{1}{8})
- (\frac{1}{2},-\frac{1}{2}, \frac{1}{2})
- \}$
-
+Note that the convection is divergence-free, $\nabla \cdot c = 0$.
<h3> Implementation </h3>
-Implementation notes:
-
-- WorkStream to parallelize local solvers. Workstream is already used in step-32, step-44.
-- Reconstructing the trace
-- Post-processing the solution for superconvergence
-- DataOutFaces: direct output of the global solution
-
+Besides implementing the above equations, the implementation below provides the following features:
+<ul>
+ <li> WorkStream to parallelize local solvers. Workstream is already used in
+ step-32, step-44.
+ <li> Reconstruct the local DG solution from the trace.
+ <li> Post-processing the solution for superconvergence.
+ <li> DataOutFaces for direct output of the global skeleton solution.
+</ul>
<h1>Results</h1>
+<h3>Program output</h3>
-When the program is run, it outputs information about the respective steps and
-convergence tables with errors in the various components in the end. In 2D,
-the convergence tables look the following:
+We first have a look at the output generated by the program when run in 2D. In
+the four images below, we show the solution for polynomial degree <i>p</i>=1
+and cycles 2, 3, 4, and 8 of the program. In the plots, we overlay the data
+generated from the internal data (DG part) with the skeleton part ($\hat{u}$)
+into the same plot. We had to generate two different data sets because cells
+and faces represent different geometric entities, the combination of which (in
+the same file) is not supported in the VTK output of deal.II.
+
+The images show the distinctive features of HDG: The cell solution (colored
+surfaces) is discontinuous between the cells. The solution on the skeleton
+variable sits on the faces and ties together the local parts. The skeleton
+solution is not continuous on the vertices where the faces meet, even though
+its values are quite close along lines in the same coordinate direction. The
+skeleton solution can be interpreted as a rubber spring between the two sides
+that balances the jumps in the solution (or rather, the flux $\kappa \nabla u
++ \mathbf{c} u$). As the mesh is refined, the jumps between the cells get
+small (we represent a smooth solution), and the skeleton solution approaches
+the interior parts. For cycle 8, there is no visible difference in the two
+variables. We also see how boundary conditions are implemented weakly and that
+the interior variables do not exactly satisfy boundary conditions. On the
+lower and left boundaries, we set Neumann boundary conditions, whereas we set
+Dirichlet conditions on the right and top boundaries.
+
+<table align="center" border="1" cellspacing="3" cellpadding="3">
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.sol_2.png" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.sol_3.png" alt="">
+ </td>
+ </tr>
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.sol_4.png" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.sol_8.png" alt="">
+ </td>
+ </tr>
+</table>
+
+Next, we have a look at the post-processed solution, again at cycles 2, 3, 4,
+and 8. This is a discontinuous solution that is locally described by second
+order polynomials. While the solution does not look very good on the mesh of
+cycle two, it looks much better for cycles three and four. As shown by the
+convergence table below, we find that is also converges more quickly to the
+analytical solution.
+
+<table align="center" border="1" cellspacing="3" cellpadding="3">
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.post_2.png" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.post_3.png" alt="">
+ </td>
+ </tr>
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.post_4.png" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.post_8.png" alt="">
+ </td>
+ </tr>
+</table>
+
+Finally, we look at the solution for <i>p</i>=3 at cycle 2. Despite the coarse
+mesh with only 64 cells, the post-processed solution is similar in quality
+to the linear solution (not post-processed) at cycle 8 with 4,096
+cells. This clearly shows the superiority of high order methods for smooth
+solutions.
+
+<table align="center" border="1" cellspacing="3" cellpadding="3">
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.sol_q3_2.png" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.post_q3_2.png" alt="">
+ </td>
+ </tr>
+</table>
+
+<h4>Convergence tables</h4>
+
+When the program is run, it also outputs information about the respective
+steps and convergence tables with errors in the various components in the
+end. In 2D, the convergence tables look the following:
@code
Q1 elements, adaptive refinement:
convergence rates of Q1 elements in the $L_2$ norm for both the scalar
variable and the gradient variable is apparent, as is the cubic rate for the
postprocessed scalar variable in the $L_2$ norm. Note that is is a distinctive
-feature of an HDG solution. In typical continuous finite element, the gradient
-of the solution of order <i>p</i> converges at rate <i>p</i> only, as opposed
-to <i>p</i>+1 for the actual solution. Even though superconvergence results
-for finite elements are also available (e.g. superconvergent patch recovery
-first introduced by Zienkiewicz and Zhu), these are typically limited to
-structured meshes and other special cases. Likewise, the scalar variable and
-gradient for Q3 elements converge at fourth order and the postprocessed scalar
+feature of an HDG solution. In typical continuous finite elements, the
+gradient of the solution of order <i>p</i> converges at rate <i>p</i> only, as
+opposed to <i>p</i>+1 for the actual solution. Even though superconvergence
+results for finite elements are also available (e.g. superconvergent patch
+recovery first introduced by Zienkiewicz and Zhu), these are typically limited
+to structured meshes and other special cases. For Q3 HDG variables, the scalar
+variable and gradient converge at fourth order and the postprocessed scalar
variable at fifth order.
The same convergence rates are observed in 3d.
110592 5419008 3.482e-05 3.94 3.055e-04 3.95 7.374e-07 5.00
@endcode
+<h3>Comparison with continuous finite elements</h3>
+
+<h4>Results for 2D</h4>
+
+The convergence tables verify the expected convergence rates stated in the
+introduction. Now, we want to show a quick comparison of the computational
+efficiency of the HDG method compared to a usual finite element (continuous
+Galkerin) method on the problem of this tutorial. Of course, stability aspects
+of the HDG method compared to continuous finite elements for
+transport-dominated problems are also important in practice, which is an
+aspect not seen on a problem with smooth analytic solution. In the picture
+below, we compare the $L_2$ error as a function of the number of degrees of
+freedom (left) and of the computing time spent in the linear solver (right)
+for two space dimensions of continuous finite elements (CG) and the hybridized
+discontinuous Galerkin method presented in this tutorial. As opposed to the
+tutorial where we only use unpreconditioned BiCGStab, the times shown in the
+figures below use the Trilinos algebraic multigrid preconditioner in
+TrilinosWrappers::PreconditionAMG. For the HDG part, a wrapper around
+ChunkSparseMatrix for the trace variable has been used in order to utilize the
+block structure in the matrix on the finest level.
+
+<table align="center" border="1" cellspacing="3" cellpadding="3">
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.2d_plain.png" width="400" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.2dt_plain.png" width="400" alt="">
+ </td>
+ </tr>
+</table>
+
+The results in the graphs show that the HDG method is slower than continuous
+finite elements at <i>p</i>=1, about equally fast for cubic elements and
+faster for sixth order elements. However, we have seen above that the HDG
+method actually produces solutions which are more accurate than what is
+represented in the original variables. Therefore, in the next two plots below
+we instead display the error of the post-processed solution for HDG (denoted
+by $p=1^*$ for example). We now see a clear advantage of HDG for the same
+amount of work for both <i>p</i>=3 and <i>p</i>=6, and about the same quality
+for <i>p</i>=1.
+
+<table align="center" border="1" cellspacing="3" cellpadding="3">
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.2d_post.png" width="400" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.2dt_post.png" width="400" alt="">
+ </td>
+ </tr>
+</table>
+
+Since the HDG method actually produces results converging as
+<i>h</i><sup><i>p</i>+2</sup>, we should compare it to a continuous Galerkin
+solution with the same asymptotic convergence behavior, i.e., FE_Q with degree
+<i>p</i>+1. If we do this, we get the convergence curves below. We see that
+CG with second order polynomials is again clearly better than HDG with
+linears. However, the advantage of HDG for higher orders remains.
+
+<table align="center" border="1" cellspacing="3" cellpadding="3">
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.2d_postb.png" width="400" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.2dt_postb.png" width="400" alt="">
+ </td>
+ </tr>
+</table>
+
+The results are in line with properties of DG methods in general: Best
+performance is typically not achieved for linear elements, but rather at
+somewhat higher order, usually around <i>p</i>=3. This is because of a
+volume-to-surface effect for discontinuous solutions with too much of the
+solution living on the surfaces and hence duplicating work when the elements
+are linear. Put in other words, DG methods are often most efficient when used
+at relatively high order, despite their focus on discontinuous (and hence,
+seemingly low accurate) representation of solutions.
+
+<h4>Resuls for 3D</h4>
+
+We now show the same figures in 3D: The first row shows the number of degrees
+of freedom and computing time versus the $L_2$ error in the scalar variable
+<i>u</i> for CG and HDG at order <i>p</i>, the second row shows the
+post-processed HDG solution instead of the original one, and the third row
+compares the post-processed HDG solution with CG at order <i>p</i>+1. In 3D,
+the volume-to-surface effect makes the cost of HDG somewhat higher and the CG
+solution is clearly better than HDG for linears by any metric. For cubics, HDG
+and CG are of similar quality, whereas HDG is again more efficient for sixth
+order polynomials. One can alternatively also use the combination of FE_DGP
+and FE_FaceP instead of (FE_DGQ, FE_FaceQ), which do not use tensor product
+polynomials of degree <i>p</i> but Legendre polynomials of <i>complete</i>
+degree <i>p</i>. There are fewer degrees of freedom on the skeleton variable
+for FE_FaceP for a given mesh size, but the solution quality (error vs. number
+of DoFs) is very similar to the results for FE_FaceQ.
+
+<table align="center" border="1" cellspacing="3" cellpadding="3">
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.3d_plain.png" width="400" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.3dt_plain.png" width="400" alt="">
+ </td>
+ </tr>
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.3d_post.png" width="400" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.3dt_post.png" width="400" alt="">
+ </td>
+ </tr>
+ <tr>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.3d_postb.png" width="400" alt="">
+ </td>
+ <td>
+ <img src="http://www.dealii.org/images/steps/developer/step-51.3dt_postb.png" width="400" alt="">
+ </td>
+ </tr>
+</table>
+
+One final note on the efficiency comparison: We tried to use general-purpose
+sparse matrix structures and similar solvers (optimal AMG preconditioners for
+both without particular tuning of the AMG parameters on any of them) to give a
+fair picture of the cost versus accuracy of two methods, on a toy example. It
+should be noted however that GMG for continuous finite elements is about a
+factor four to five faster for <i>p</i>=3 and <i>p</i>=6. The authors of this
+tutorial have not seen similarly advanced solvers for the HDG linear
+systems. Also, there are other implementation aspects for CG available such as
+fast matrix-free approaches as shown in step-37 that make higher order
+continuous elements more competitive. Again, it is not clear to the authors of
+the tutorial whether similar improvements could be made for HDG.
<h3>Possibilities for improvements</h3>
-This tutorial program shows the a basic implementation of an HDG method
-without particular focus on a high performance implementation. Let us look at
+As already mentioned in the introduction, one possibility is to implement
+another post-processing technique as discussed in the literature.
+
+A second item that is not done optimally relates to the performance of this
+program, which is of course an issue in practical applications (weighing in
+also the better solution quality of (H)DG methods for transport-dominated
+problems). Let us look at
the computing time of the tutorial program and the share of the individual
components:
multigrid preconditioner from Trilinos. For diffusion-dominated problems as
the problem at hand with finer meshes, such a solver can be designed that
uses the matrix-vector products from the more efficient ChunkSparseMatrix on
- the finest level, as long as we are not working in parallel with MPI.
+ the finest level, as long as we are not working in parallel with MPI. For
+ MPI-parallelized computation, a standard TrilinosWrappers::SparseMatrix can
+ be used.
<li> Speed up assembly by pre-assembling parts that do not change from one
cell to another (those that do neither contain variable coefficients nor
// the simulation.
#include <deal.II/numerics/data_out_faces.h>
+
+// We start by putting the class into its own namespace.
namespace Step51
{
-using namespace dealii;
+ using namespace dealii;
// @sect3{Equation data}
//
// two exceptions. Firstly, we also create a solution for the 3d case, and
// secondly, we scale the solution so its norm is of order unity for all
// values of the solution width.
-template <int dim>
-class SolutionBase
-{
-protected:
- static const unsigned int n_source_centers = 3;
- static const Point<dim> source_centers[n_source_centers];
- static const double width;
-};
-
-
-template <>
-const Point<1>
-SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
-= { Point<1>(-1.0 / 3.0),
- Point<1>(0.0),
- Point<1>(+1.0 / 3.0)
+ template <int dim>
+ class SolutionBase
+ {
+ protected:
+ static const unsigned int n_source_centers = 3;
+ static const Point<dim> source_centers[n_source_centers];
+ static const double width;
};
-template <>
-const Point<2>
-SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
-= { Point<2>(-0.5, +0.5),
- Point<2>(-0.5, -0.5),
- Point<2>(+0.5, -0.5)
- };
+ template <>
+ const Point<1>
+ SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
+ = { Point<1>(-1.0 / 3.0),
+ Point<1>(0.0),
+ Point<1>(+1.0 / 3.0)
+ };
-template <>
-const Point<3>
-SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers]
-= { Point<3>(-0.5, +0.5, 0.25),
- Point<3>(-0.6, -0.5, -0.125),
- Point<3>(+0.5, -0.5, 0.5)
- };
-template <int dim>
-const double SolutionBase<dim>::width = 1./5.;
+ template <>
+ const Point<2>
+ SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
+ = { Point<2>(-0.5, +0.5),
+ Point<2>(-0.5, -0.5),
+ Point<2>(+0.5, -0.5)
+ };
+ template <>
+ const Point<3>
+ SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers]
+ = { Point<3>(-0.5, +0.5, 0.25),
+ Point<3>(-0.6, -0.5, -0.125),
+ Point<3>(+0.5, -0.5, 0.5)
+ };
-template <int dim>
-class Solution : public Function<dim>,
- protected SolutionBase<dim>
-{
-public:
- Solution () : Function<dim>() {}
+ template <int dim>
+ const double SolutionBase<dim>::width = 1./5.;
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
- virtual Tensor<1,dim> gradient (const Point<dim> &p,
- const unsigned int component = 0) const;
-};
+ template <int dim>
+ class Solution : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ Solution () : Function<dim>() {}
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+ };
-template <int dim>
-double Solution<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- double return_value = 0;
- for (unsigned int i=0; i<this->n_source_centers; ++i)
- {
- const Point<dim> x_minus_xi = p - this->source_centers[i];
- return_value += std::exp(-x_minus_xi.square() /
- (this->width * this->width));
- }
- return return_value /
- Utilities::fixed_power<dim>(std::sqrt(2. * numbers::PI) * this->width);
-}
+ template <int dim>
+ double Solution<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ double return_value = 0;
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> x_minus_xi = p - this->source_centers[i];
+ return_value += std::exp(-x_minus_xi.square() /
+ (this->width * this->width));
+ }
+ return return_value /
+ Utilities::fixed_power<dim>(std::sqrt(2. * numbers::PI) * this->width);
+ }
-template <int dim>
-Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &p,
- const unsigned int) const
-{
- Tensor<1,dim> return_value;
- for (unsigned int i=0; i<this->n_source_centers; ++i)
- {
- const Point<dim> x_minus_xi = p - this->source_centers[i];
- return_value += (-2 / (this->width * this->width) *
- std::exp(-x_minus_xi.square() /
- (this->width * this->width)) *
- x_minus_xi);
- }
+ template <int dim>
+ Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &p,
+ const unsigned int) const
+ {
+ Tensor<1,dim> return_value;
- return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI) *
- this->width);
-}
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> x_minus_xi = p - this->source_centers[i];
+
+ return_value += (-2 / (this->width * this->width) *
+ std::exp(-x_minus_xi.square() /
+ (this->width * this->width)) *
+ x_minus_xi);
+ }
+
+ return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI) *
+ this->width);
+ }
// gradient are collected together. This function is used when computing the
// error of the HDG approximation and its implementation is to simply call
// value and gradient function of the Solution class.
-template <int dim>
-class SolutionAndGradient : public Function<dim>,
- protected SolutionBase<dim>
-{
-public:
- SolutionAndGradient () : Function<dim>(dim) {}
+ template <int dim>
+ class SolutionAndGradient : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ SolutionAndGradient () : Function<dim>(dim) {}
- virtual void vector_value (const Point<dim> &p,
- Vector<double> &v) const;
-};
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const;
+ };
-template <int dim>
-void SolutionAndGradient<dim>::vector_value (const Point<dim> &p,
- Vector<double> &v) const
-{
- AssertDimension(v.size(), dim+1);
- Solution<dim> solution;
- Tensor<1,dim> grad = solution.gradient(p);
- for (unsigned int d=0; d<dim; ++d)
- v[d] = -grad[d];
- v[dim] = solution.value(p);
-}
+ template <int dim>
+ void SolutionAndGradient<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ AssertDimension(v.size(), dim+1);
+ Solution<dim> solution;
+ Tensor<1,dim> grad = solution.gradient(p);
+ for (unsigned int d=0; d<dim; ++d)
+ v[d] = -grad[d];
+ v[dim] = solution.value(p);
+ }
// Next comes the implementation of the convection velocity. As described in
// the introduction, we choose a velocity field that is $(y, -x)$ in 2D and
// $(y, -x, 1)$ in 3D. This gives a divergence-free velocity field.
-template <int dim>
-class ConvectionVelocity : public TensorFunction<1,dim>
-{
-public:
- ConvectionVelocity() : TensorFunction<1,dim>() {}
+ template <int dim>
+ class ConvectionVelocity : public TensorFunction<1,dim>
+ {
+ public:
+ ConvectionVelocity() : TensorFunction<1,dim>() {}
- virtual Tensor<1,dim> value (const Point<dim> &p) const;
-};
+ virtual Tensor<1,dim> value (const Point<dim> &p) const;
+ };
-template <int dim>
-Tensor<1,dim>
-ConvectionVelocity<dim>::value(const Point<dim> &p) const
-{
- Tensor<1,dim> convection;
- switch (dim)
- {
- case 1:
- convection[0] = 1;
- break;
- case 2:
- convection[0] = p[1];
- convection[1] = -p[0];
- break;
- case 3:
- convection[0] = p[1];
- convection[1] = -p[0];
- convection[2] = 1;
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
- return convection;
-}
+ template <int dim>
+ Tensor<1,dim>
+ ConvectionVelocity<dim>::value(const Point<dim> &p) const
+ {
+ Tensor<1,dim> convection;
+ switch (dim)
+ {
+ case 1:
+ convection[0] = 1;
+ break;
+ case 2:
+ convection[0] = p[1];
+ convection[1] = -p[0];
+ break;
+ case 3:
+ convection[0] = p[1];
+ convection[1] = -p[0];
+ convection[2] = 1;
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ return convection;
+ }
// a convection term instead of the reaction term. Since the velocity field is
// incompressible, i.e. $\nabla \cdot \mathbf{c} = 0$, this term simply reads
// $\mathbf{c} \nabla \ve u$.
-template <int dim>
-class RightHandSide : public Function<dim>,
- protected SolutionBase<dim>
-{
-public:
- RightHandSide () : Function<dim>() {}
+ template <int dim>
+ class RightHandSide : public Function<dim>,
+ protected SolutionBase<dim>
+ {
+ public:
+ RightHandSide () : Function<dim>() {}
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
-private:
- const ConvectionVelocity<dim> convection_velocity;
-};
+ private:
+ const ConvectionVelocity<dim> convection_velocity;
+ };
-template <int dim>
-double RightHandSide<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- Tensor<1,dim> convection = convection_velocity.value(p);
- double return_value = 0;
- for (unsigned int i=0; i<this->n_source_centers; ++i)
- {
- const Point<dim> x_minus_xi = p - this->source_centers[i];
-
- return_value +=
- ((2*dim - 2*convection*x_minus_xi - 4*x_minus_xi.square()/
- (this->width * this->width)) /
- (this->width * this->width) *
- std::exp(-x_minus_xi.square() /
- (this->width * this->width)));
- }
+ template <int dim>
+ double RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+ {
+ Tensor<1,dim> convection = convection_velocity.value(p);
+ double return_value = 0;
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> x_minus_xi = p - this->source_centers[i];
+
+ return_value +=
+ ((2*dim - 2*convection*x_minus_xi - 4*x_minus_xi.square()/
+ (this->width * this->width)) /
+ (this->width * this->width) *
+ std::exp(-x_minus_xi.square() /
+ (this->width * this->width)));
+ }
- return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI)
- * this->width);
-}
+ return return_value / Utilities::fixed_power<dim>(std::sqrt(2 * numbers::PI)
+ * this->width);
+ }
// @sect3{The HDG solver class}
// generate the system matrix and once when we compute the element-interior
// solutions from the skeleton values) and for the postprocessing where
// we extract a solution that converges at higher order.
-template <int dim>
-class HDG
-{
-public:
- enum RefinementMode
+ template <int dim>
+ class HDG
{
- global_refinement, adaptive_refinement
+ public:
+ enum RefinementMode
+ {
+ global_refinement, adaptive_refinement
+ };
+
+ HDG (const unsigned int degree,
+ const RefinementMode refinement_mode);
+ void run ();
+
+ private:
+
+ void setup_system ();
+ void assemble_system (const bool reconstruct_trace = false);
+ void solve ();
+ void postprocess ();
+ void refine_grid (const unsigned int cylce);
+ void output_results (const unsigned int cycle);
+
+ // Data for the assembly and solution of the primal variables.
+ struct PerTaskData;
+ struct ScratchData;
+
+ // Post-processing the solution to obtain $u^*$ is an element-by-element
+ // procedure; as such, we do not need to assemble any global data and do
+ // not declare any 'task data' for WorkStream to use.
+ struct PostProcessScratchData;
+
+ // The following three functions are used by WorkStream to do the actual
+ // work of the program.
+ void assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ ScratchData &scratch,
+ PerTaskData &task_data);
+
+ void copy_local_to_global(const PerTaskData &data);
+
+ void postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ PostProcessScratchData &scratch,
+ unsigned int &empty_data);
+
+
+ Triangulation<dim> triangulation;
+
+ // The 'local' solutions are interior to each element. These
+ // represent the primal solution field $u$ as well as the auxiliary
+ // field $\mathbf{q}$.
+ FESystem<dim> fe_local;
+ DoFHandler<dim> dof_handler_local;
+ Vector<double> solution_local;
+
+ // The new finite element type and corresponding <code>DoFHandler</code> are
+ // used for the global skeleton solution that couples the element-level local
+ // solutions.
+ FE_FaceQ<dim> fe;
+ DoFHandler<dim> dof_handler;
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ // As stated in the introduction, HDG solutions can be post-processed to
+ // attain superconvergence rates of $\mathcal{O}(h^{p+2})$. The
+ // post-processed solution is a discontinuous finite element solution
+ // representing the primal variable on the interior of each cell. We define
+ // a FE type of degree $p+1$ to represent this post-processed solution,
+ // which we only use for output after constructing it.
+ FE_DGQ<dim> fe_u_post;
+ DoFHandler<dim> dof_handler_u_post;
+ Vector<double> solution_u_post;
+
+ // The degrees of freedom corresponding to the skeleton strongly enforce
+ // Dirichlet boundary conditions, just as in a continuous Galerkin finite
+ // element method. We can enforce the boundary conditions in an analogous
+ // manner through the use of <code>ConstrainMatrix</code> constructs. In
+ // addition, hanging nodes are handled in the same way as for
+ // continuous finite elements: For the face elements which
+ // only define degrees of freedom on the face, this process sets the
+ // solution on the refined to be the one from the coarse side.
+ ConstraintMatrix constraints;
+
+ // The usage of the ChunkSparseMatrix class is similar to the usual sparse
+ // matrices: You need a sparsity pattern of type ChunkSparsityPattern and
+ // the actual matrix object. When creating the sparsity pattern, we just
+ // have to additionally pass the size of local blocks.
+ ChunkSparsityPattern sparsity_pattern;
+ ChunkSparseMatrix<double> system_matrix;
+
+ // Same as step-7:
+ const RefinementMode refinement_mode;
+ ConvergenceTable convergence_table;
};
- HDG (const unsigned int degree,
- const RefinementMode refinement_mode);
- void run ();
-
-private:
-
-// Data for the assembly and solution of the primal variables.
- struct PerTaskData;
- struct ScratchData;
-
-// Post-processing the solution to obtain $u^*$ is an element-by-element
-// procedure; as such, we do not need to assemble any global data and do
-// not declare any 'task data' for WorkStream to use.
- struct PostProcessScratchData;
-
- void setup_system ();
- void assemble_system (const bool reconstruct_trace = false);
- void solve ();
- void postprocess ();
-
- void refine_grid (const unsigned int cylce);
- void output_results (const unsigned int cycle);
-
-// The following three functions are used by WorkStream to do the actual work of
-// the program.
- void assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
- ScratchData &scratch,
- PerTaskData &task_data);
-
- void copy_local_to_global(const PerTaskData &data);
-
- void postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
- PostProcessScratchData &scratch,
- unsigned int &empty_data);
-
-
- Triangulation<dim> triangulation;
-
- // The 'local' solutions are interior to each element. These
- // represent the primal solution field $u$ as well as the auxiliary
- // field $\mathbf{q} = -\nabla u$.
- FESystem<dim> fe_local;
- DoFHandler<dim> dof_handler_local;
- Vector<double> solution_local;
-
- // The new finite element type and corresponding <code>DoFHandler</code> are
- // used for the global skeleton solution that couples the element-level local
- // solutions.
- FE_FaceQ<dim> fe;
- DoFHandler<dim> dof_handler;
- Vector<double> solution;
- Vector<double> system_rhs;
-
- // As stated in the introduction, HDG solutions can be post-processed to
- // attain superconvergence rates of $\mathcal{O}(h^{p+2})$. The
- // post-processed solution is a discontinuous finite element solution
- // representing the primal variable on the interior of each cell. We define
- // a FE type of degree $p+1$ to represent this post-processed solution,
- // which we only use for output after constructing it.
- FE_DGQ<dim> fe_u_post;
- DoFHandler<dim> dof_handler_u_post;
- Vector<double> solution_u_post;
-
- // The degrees of freedom corresponding to the skeleton strongly enforce
- // Dirichlet boundary conditions, just as in a continuous Galerkin finite
- // element method. We can enforce the boundary conditions in an analogous
- // manner through the use of <code>ConstrainMatrix</code> constructs. In
- // addition, hanging nodes where cells of different refinement levels meet
- // are set as for continuous finite elements: For the face elements which
- // only define degrees of freedom on the face, this process sets the
- // solution on the refined to be the one from the coarse side.
- ConstraintMatrix constraints;
-
- // The usage of the ChunkSparseMatrix class is similar to the usual sparse
- // matrices: You need a sparsity pattern of type ChunkSparsityPattern and
- // the actual matrix object. When creating the sparsity pattern, we just
- // have to additionally pass the size of local blocks.
- ChunkSparsityPattern sparsity_pattern;
- ChunkSparseMatrix<double> system_matrix;
-
- // Same as step-7:
- const RefinementMode refinement_mode;
- ConvergenceTable convergence_table;
-};
-
-// @sect3{The HDG class implementation}
-
-// @sect4{Constructor} The constructor is similar to those in other examples,
-// with the exception of handling multiple <code>DoFHandler</code> and
-// <code>FiniteElement</code> objects. Note that we create a system of finite
-// elements for the local DG part, including the gradient/flux part and the
-// scalar part.
-template <int dim>
-HDG<dim>::HDG (const unsigned int degree,
- const RefinementMode refinement_mode) :
- fe_local (FE_DGQ<dim>(degree), dim,
- FE_DGQ<dim>(degree), 1),
- dof_handler_local (triangulation),
- fe (degree),
- dof_handler (triangulation),
- fe_u_post (degree+1),
- dof_handler_u_post (triangulation),
- refinement_mode (refinement_mode)
-{}
-
-
-
-// @sect4{HDG::setup_system}
-// The system for an HDG solution is setup in an analogous manner to most
-// of the other tutorial programs. We are careful to distribute dofs with
-// all of our <code>DoFHandler</code> objects. The @p solution and @p system_matrix
-// objects go with the global skeleton solution.
-template <int dim>
-void
-HDG<dim>::setup_system ()
-{
- dof_handler_local.distribute_dofs(fe_local);
- dof_handler.distribute_dofs(fe);
- dof_handler_u_post.distribute_dofs(fe_u_post);
-
- std::cout << " Number of degrees of freedom: "
- << dof_handler.n_dofs()
- << std::endl;
-
- solution.reinit (dof_handler.n_dofs());
- system_rhs.reinit (dof_handler.n_dofs());
-
- solution_local.reinit (dof_handler_local.n_dofs());
- solution_u_post.reinit (dof_handler_u_post.n_dofs());
-
- constraints.clear ();
- DoFTools::make_hanging_node_constraints (dof_handler, constraints);
- typename FunctionMap<dim>::type boundary_functions;
- Solution<dim> solution_function;
- boundary_functions[0] = &solution_function;
- VectorTools::project_boundary_values (dof_handler,
- boundary_functions,
- QGauss<dim-1>(fe.degree+1),
- constraints);
- constraints.close ();
+ // @sect3{The HDG class implementation}
+
+ // @sect4{Constructor}
+ // The constructor is similar to those in other examples,
+ // with the exception of handling multiple <code>DoFHandler</code> and
+ // <code>FiniteElement</code> objects. Note that we create a system of finite
+ // elements for the local DG part, including the gradient/flux part and the
+ // scalar part.
+ template <int dim>
+ HDG<dim>::HDG (const unsigned int degree,
+ const RefinementMode refinement_mode) :
+ fe_local (FE_DGQ<dim>(degree), dim,
+ FE_DGQ<dim>(degree), 1),
+ dof_handler_local (triangulation),
+ fe (degree),
+ dof_handler (triangulation),
+ fe_u_post (degree+1),
+ dof_handler_u_post (triangulation),
+ refinement_mode (refinement_mode)
+ {}
- {
- CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern (dof_handler, csp,
- constraints, false);
- sparsity_pattern.copy_from(csp, fe.dofs_per_face);
- }
- system_matrix.reinit (sparsity_pattern);
-}
+ // @sect4{HDG::setup_system}
+ // The system for an HDG solution is setup in an analogous manner to most
+ // of the other tutorial programs. We are careful to distribute dofs with
+ // all of our <code>DoFHandler</code> objects. The @p solution and @p system_matrix
+ // objects go with the global skeleton solution.
+ template <int dim>
+ void
+ HDG<dim>::setup_system ()
+ {
+ dof_handler_local.distribute_dofs(fe_local);
+ dof_handler.distribute_dofs(fe);
+ dof_handler_u_post.distribute_dofs(fe_u_post);
+
+ std::cout << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+
+ solution_local.reinit (dof_handler_local.n_dofs());
+ solution_u_post.reinit (dof_handler_u_post.n_dofs());
+
+ constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+ typename FunctionMap<dim>::type boundary_functions;
+ Solution<dim> solution_function;
+ boundary_functions[0] = &solution_function;
+ VectorTools::project_boundary_values (dof_handler,
+ boundary_functions,
+ QGauss<dim-1>(fe.degree+1),
+ constraints);
+ constraints.close ();
+
+ // When creating the chunk sparsity pattern, we first create the usual
+ // compressed sparsity pattern and then set the chunk size, which is equal
+ // to the number of dofs on a face, when copying this into the final
+ // sparsity pattern.
+ {
+ CompressedSimpleSparsityPattern csp (dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, csp,
+ constraints, false);
+ sparsity_pattern.copy_from(csp, fe.dofs_per_face);
+ }
+ system_matrix.reinit (sparsity_pattern);
+ }
-// @sect4{HDG::PerTaskData} Next come the definition of the local data
-// structures for the parallel assembly. The first structure @p PerTaskData
-// contains the local vector and matrix that are written into the global
-// matrix, whereas the ScratchData contains all data that we need for the
-// local assembly. There is one variable worth noting here, namely the boolean
-// variable @p trace_reconstruct. As mentioned introdution, we solve the HDG
-// system in two steps. First, we create a linear system for the skeleton
-// system where we condense the local part into it by $D-CA^{-1}B$. Then, we
-// solve for the local part using the skeleton solution. For these two steps,
-// we need the same matrices on the elements twice, which we want to compute
-// by two assembly steps. Since most of the code is similar, we do this with
-// the same function but only switch between the two based on a flag that we
-// set when starting the assembly. Since we need to pass this information on
-// to the local worker routines, we store it once in the task data.
-template <int dim>
-struct HDG<dim>::PerTaskData
-{
- FullMatrix<double> cell_matrix;
- Vector<double> cell_vector;
- std::vector<types::global_dof_index> dof_indices;
- bool trace_reconstruct;
- PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct)
- : cell_matrix(n_dofs, n_dofs),
- cell_vector(n_dofs),
- dof_indices(n_dofs),
- trace_reconstruct(trace_reconstruct)
- {}
-};
-
-
-
-// @sect4{HDG::ScratchData} @p ScratchData contains persistent data for each
-// thread within <code>WorkStream</code>. The <code>FEValues</code>, matrix,
-// and vector objects should be familiar by now. There are two objects that
-// need to be discussed: @p std::vector<std::vector<unsigned int> >
-// fe_local_support_on_face and @p std::vector<std::vector<unsigned int> >
-// fe_support_on_face. These are used to indicate whether or not the finite
-// elements chosen have support (non-zero values) on a given face of the
-// reference cell for the local part associated to @p fe_local and the
-// skeleton part @p f, which is why we can extract this information in the
-// constructor and store it once for all cells that we work on. Had we not
-// stored this information, we would be forced to assemble a large number of
-// zero terms on each cell, which would significantly slow the program.
-template <int dim>
-struct HDG<dim>::ScratchData
-{
- FEValues<dim> fe_values_local;
- FEFaceValues<dim> fe_face_values_local;
- FEFaceValues<dim> fe_face_values;
-
- FullMatrix<double> ll_matrix;
- FullMatrix<double> lf_matrix;
- FullMatrix<double> fl_matrix;
- FullMatrix<double> tmp_matrix;
- Vector<double> l_rhs;
- Vector<double> tmp_rhs;
-
- std::vector<Tensor<1,dim> > q_phi;
- std::vector<double> q_phi_div;
- std::vector<double> u_phi;
- std::vector<Tensor<1,dim> > u_phi_grad;
- std::vector<double> tr_phi;
- std::vector<double> trace_values;
-
- std::vector<std::vector<unsigned int> > fe_local_support_on_face;
- std::vector<std::vector<unsigned int> > fe_support_on_face;
-
- ConvectionVelocity<dim> convection_velocity;
- RightHandSide<dim> right_hand_side;
- const Solution<dim> exact_solution;
-
- ScratchData(const FiniteElement<dim> &fe,
- const FiniteElement<dim> &fe_local,
- const QGauss<dim> &quadrature_formula,
- const QGauss<dim-1> &face_quadrature_formula,
- const UpdateFlags local_flags,
- const UpdateFlags local_face_flags,
- const UpdateFlags flags)
- :
- fe_values_local (fe_local, quadrature_formula, local_flags),
- fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags),
- fe_face_values (fe, face_quadrature_formula, flags),
- ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell),
- lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell),
- fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
- tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
- l_rhs (fe_local.dofs_per_cell),
- tmp_rhs (fe_local.dofs_per_cell),
- q_phi (fe_local.dofs_per_cell),
- q_phi_div (fe_local.dofs_per_cell),
- u_phi (fe_local.dofs_per_cell),
- u_phi_grad (fe_local.dofs_per_cell),
- tr_phi (fe.dofs_per_cell),
- trace_values(face_quadrature_formula.size()),
- fe_local_support_on_face(GeometryInfo<dim>::faces_per_cell),
- fe_support_on_face(GeometryInfo<dim>::faces_per_cell)
+ // @sect4{HDG::PerTaskData}
+ // Next comes the definition of the local data structures for the parallel
+ // assembly. The first structure @p PerTaskData contains the local vector
+ // and matrix that are written into the global matrix, whereas the
+ // ScratchData contains all data that we need for the local assembly. There
+ // is one variable worth noting here, namely the boolean variable @p
+ // trace_reconstruct. As mentioned in the introdution, we solve the HDG
+ // system in two steps. First, we create a linear system for the skeleton
+ // system where we condense the local part into it via the Schur complement
+ // $D-CA^{-1}B$. Then, we solve for the local part using the skeleton
+ // solution. For these two steps, we need the same matrices on the elements
+ // twice, which we want to compute by two assembly steps. Since most of the
+ // code is similar, we do this with the same function but only switch
+ // between the two based on a flag that we set when starting the
+ // assembly. Since we need to pass this information on to the local worker
+ // routines, we store it once in the task data.
+ template <int dim>
+ struct HDG<dim>::PerTaskData
{
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- for (unsigned int i=0; i<fe_local.dofs_per_cell; ++i)
- {
- if (fe_local.has_support_on_face(i,face))
- fe_local_support_on_face[face].push_back(i);
- }
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_vector;
+ std::vector<types::global_dof_index> dof_indices;
+
+ bool trace_reconstruct;
+
+ PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct)
+ : cell_matrix(n_dofs, n_dofs),
+ cell_vector(n_dofs),
+ dof_indices(n_dofs),
+ trace_reconstruct(trace_reconstruct)
+ {}
+ };
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
- {
- if (fe.has_support_on_face(i,face))
- fe_support_on_face[face].push_back(i);
- }
- }
- ScratchData(const ScratchData &sd)
- :
- fe_values_local (sd.fe_values_local.get_fe(),
- sd.fe_values_local.get_quadrature(),
- sd.fe_values_local.get_update_flags()),
- fe_face_values_local (sd.fe_face_values_local.get_fe(),
- sd.fe_face_values_local.get_quadrature(),
- sd.fe_face_values_local.get_update_flags()),
- fe_face_values (sd.fe_face_values.get_fe(),
- sd.fe_face_values.get_quadrature(),
- sd.fe_face_values.get_update_flags()),
- ll_matrix (sd.ll_matrix),
- lf_matrix (sd.lf_matrix),
- fl_matrix (sd.fl_matrix),
- tmp_matrix (sd.tmp_matrix),
- l_rhs (sd.l_rhs),
- tmp_rhs (sd.tmp_rhs),
- q_phi (sd.q_phi),
- q_phi_div (sd.q_phi_div),
- u_phi (sd.u_phi),
- u_phi_grad (sd.u_phi_grad),
- tr_phi (sd.tr_phi),
- trace_values(sd.trace_values),
- fe_local_support_on_face(sd.fe_local_support_on_face),
- fe_support_on_face(sd.fe_support_on_face)
- {}
-};
+ // @sect4{HDG::ScratchData}
+ // @p ScratchData contains persistent data for each
+ // thread within <code>WorkStream</code>. The <code>FEValues</code>, matrix,
+ // and vector objects should be familiar by now. There are two objects that
+ // need to be discussed: @p std::vector<std::vector<unsigned int> >
+ // fe_local_support_on_face and @p std::vector<std::vector<unsigned int> >
+ // fe_support_on_face. These are used to indicate whether or not the finite
+ // elements chosen have support (non-zero values) on a given face of the
+ // reference cell for the local part associated to @p fe_local and the
+ // skeleton part @p fe. We extract this information in the
+ // constructor and store it once for all cells that we work on. Had we not
+ // stored this information, we would be forced to assemble a large number of
+ // zero terms on each cell, which would significantly slow the program.
+ template <int dim>
+ struct HDG<dim>::ScratchData
+ {
+ FEValues<dim> fe_values_local;
+ FEFaceValues<dim> fe_face_values_local;
+ FEFaceValues<dim> fe_face_values;
+
+ FullMatrix<double> ll_matrix;
+ FullMatrix<double> lf_matrix;
+ FullMatrix<double> fl_matrix;
+ FullMatrix<double> tmp_matrix;
+ Vector<double> l_rhs;
+ Vector<double> tmp_rhs;
+
+ std::vector<Tensor<1,dim> > q_phi;
+ std::vector<double> q_phi_div;
+ std::vector<double> u_phi;
+ std::vector<Tensor<1,dim> > u_phi_grad;
+ std::vector<double> tr_phi;
+ std::vector<double> trace_values;
+
+ std::vector<std::vector<unsigned int> > fe_local_support_on_face;
+ std::vector<std::vector<unsigned int> > fe_support_on_face;
+
+ ConvectionVelocity<dim> convection_velocity;
+ RightHandSide<dim> right_hand_side;
+ const Solution<dim> exact_solution;
+
+ ScratchData(const FiniteElement<dim> &fe,
+ const FiniteElement<dim> &fe_local,
+ const QGauss<dim> &quadrature_formula,
+ const QGauss<dim-1> &face_quadrature_formula,
+ const UpdateFlags local_flags,
+ const UpdateFlags local_face_flags,
+ const UpdateFlags flags)
+ :
+ fe_values_local (fe_local, quadrature_formula, local_flags),
+ fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags),
+ fe_face_values (fe, face_quadrature_formula, flags),
+ ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell),
+ lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell),
+ fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
+ tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell),
+ l_rhs (fe_local.dofs_per_cell),
+ tmp_rhs (fe_local.dofs_per_cell),
+ q_phi (fe_local.dofs_per_cell),
+ q_phi_div (fe_local.dofs_per_cell),
+ u_phi (fe_local.dofs_per_cell),
+ u_phi_grad (fe_local.dofs_per_cell),
+ tr_phi (fe.dofs_per_cell),
+ trace_values(face_quadrature_formula.size()),
+ fe_local_support_on_face(GeometryInfo<dim>::faces_per_cell),
+ fe_support_on_face(GeometryInfo<dim>::faces_per_cell)
+ {
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int i=0; i<fe_local.dofs_per_cell; ++i)
+ {
+ if (fe_local.has_support_on_face(i,face))
+ fe_local_support_on_face[face].push_back(i);
+ }
+
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ {
+ if (fe.has_support_on_face(i,face))
+ fe_support_on_face[face].push_back(i);
+ }
+ }
+ ScratchData(const ScratchData &sd)
+ :
+ fe_values_local (sd.fe_values_local.get_fe(),
+ sd.fe_values_local.get_quadrature(),
+ sd.fe_values_local.get_update_flags()),
+ fe_face_values_local (sd.fe_face_values_local.get_fe(),
+ sd.fe_face_values_local.get_quadrature(),
+ sd.fe_face_values_local.get_update_flags()),
+ fe_face_values (sd.fe_face_values.get_fe(),
+ sd.fe_face_values.get_quadrature(),
+ sd.fe_face_values.get_update_flags()),
+ ll_matrix (sd.ll_matrix),
+ lf_matrix (sd.lf_matrix),
+ fl_matrix (sd.fl_matrix),
+ tmp_matrix (sd.tmp_matrix),
+ l_rhs (sd.l_rhs),
+ tmp_rhs (sd.tmp_rhs),
+ q_phi (sd.q_phi),
+ q_phi_div (sd.q_phi_div),
+ u_phi (sd.u_phi),
+ u_phi_grad (sd.u_phi_grad),
+ tr_phi (sd.tr_phi),
+ trace_values(sd.trace_values),
+ fe_local_support_on_face(sd.fe_local_support_on_face),
+ fe_support_on_face(sd.fe_support_on_face)
+ {}
+ };
-// @sect4{HDG::PostProcessScratchData}
-// @p PostProcessScratchData contains the data used by <code>WorkStream</code>
-// when post-processing the local solution $u^*$. It is similar, but much
-// simpler, than @p ScratchData.
-template <int dim>
-struct HDG<dim>::PostProcessScratchData
-{
- FEValues<dim> fe_values_local;
- FEValues<dim> fe_values;
-
- std::vector<double> u_values;
- std::vector<Tensor<1,dim> > u_gradients;
- FullMatrix<double> cell_matrix;
-
- Vector<double> cell_rhs;
- Vector<double> cell_sol;
-
- PostProcessScratchData(const FiniteElement<dim> &fe,
- const FiniteElement<dim> &fe_local,
- const QGauss<dim> &quadrature_formula,
- const UpdateFlags local_flags,
- const UpdateFlags flags)
- :
- fe_values_local (fe_local, quadrature_formula, local_flags),
- fe_values (fe, quadrature_formula, flags),
- u_values (quadrature_formula.size()),
- u_gradients (quadrature_formula.size()),
- cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell),
- cell_rhs (fe.dofs_per_cell),
- cell_sol (fe.dofs_per_cell)
- {}
- PostProcessScratchData(const PostProcessScratchData &sd)
- :
- fe_values_local (sd.fe_values_local.get_fe(),
- sd.fe_values_local.get_quadrature(),
- sd.fe_values_local.get_update_flags()),
- fe_values (sd.fe_values.get_fe(),
- sd.fe_values.get_quadrature(),
- sd.fe_values.get_update_flags()),
- u_values (sd.u_values),
- u_gradients (sd.u_gradients),
- cell_matrix (sd.cell_matrix),
- cell_rhs (sd.cell_rhs),
- cell_sol (sd.cell_sol)
- {}
-};
+ // @sect4{HDG::PostProcessScratchData}
+ // @p PostProcessScratchData contains the data used by <code>WorkStream</code>
+ // when post-processing the local solution $u^*$. It is similar, but much
+ // simpler, than @p ScratchData.
+ template <int dim>
+ struct HDG<dim>::PostProcessScratchData
+ {
+ FEValues<dim> fe_values_local;
+ FEValues<dim> fe_values;
+
+ std::vector<double> u_values;
+ std::vector<Tensor<1,dim> > u_gradients;
+ FullMatrix<double> cell_matrix;
+
+ Vector<double> cell_rhs;
+ Vector<double> cell_sol;
+
+ PostProcessScratchData(const FiniteElement<dim> &fe,
+ const FiniteElement<dim> &fe_local,
+ const QGauss<dim> &quadrature_formula,
+ const UpdateFlags local_flags,
+ const UpdateFlags flags)
+ :
+ fe_values_local (fe_local, quadrature_formula, local_flags),
+ fe_values (fe, quadrature_formula, flags),
+ u_values (quadrature_formula.size()),
+ u_gradients (quadrature_formula.size()),
+ cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell),
+ cell_rhs (fe.dofs_per_cell),
+ cell_sol (fe.dofs_per_cell)
+ {}
+
+ PostProcessScratchData(const PostProcessScratchData &sd)
+ :
+ fe_values_local (sd.fe_values_local.get_fe(),
+ sd.fe_values_local.get_quadrature(),
+ sd.fe_values_local.get_update_flags()),
+ fe_values (sd.fe_values.get_fe(),
+ sd.fe_values.get_quadrature(),
+ sd.fe_values.get_update_flags()),
+ u_values (sd.u_values),
+ u_gradients (sd.u_gradients),
+ cell_matrix (sd.cell_matrix),
+ cell_rhs (sd.cell_rhs),
+ cell_sol (sd.cell_sol)
+ {}
+ };
-// @sect4{HDG::copy_local_to_global}
-// If we are in the first step of the solution, i.e. @p trace_reconstruct=false,
-// then we assemble the global system.
-template <int dim>
-void HDG<dim>::copy_local_to_global(const PerTaskData &data)
-{
- if (data.trace_reconstruct == false)
- constraints.distribute_local_to_global (data.cell_matrix,
- data.cell_vector,
- data.dof_indices,
- system_matrix, system_rhs);
-}
+ // @sect4{HDG::assemble_system}
+ // The @p assemble_system function is similar to <code>Step-32</code>, where
+ // the quadrature formula and the update flags are set up, and then
+ // <code>WorkStream</code> is used to do the work in a multi-threaded
+ // manner. The @p trace_reconstruct input parameter is used to decide
+ // whether we are solving for the global skeleton solution (false) or the
+ // local solution (true).
+ template <int dim>
+ void
+ HDG<dim>::assemble_system (const bool trace_reconstruct)
+ {
+ const QGauss<dim> quadrature_formula(fe.degree+1);
+ const QGauss<dim-1> face_quadrature_formula(fe.degree+1);
+ const UpdateFlags local_flags (update_values | update_gradients |
+ update_JxW_values | update_quadrature_points);
-// @sect4{HDG::assemble_system}
-// The @p assemble_system function is similar to <code>Step-32</code>, where
-// the quadrature formula and the update flags are set up, and then
-// <code>WorkStream</code> is used to do the work in a multi-threaded manner.
-// The @p trace_reconstruct input parameter is used to decide whether we are
-// solving for the local solution (true) or the global skeleton solution
-// (false).
-template <int dim>
-void
-HDG<dim>::assemble_system (const bool trace_reconstruct)
-{
- const QGauss<dim> quadrature_formula(fe.degree+1);
- const QGauss<dim-1> face_quadrature_formula(fe.degree+1);
-
- const UpdateFlags local_flags (update_values | update_gradients |
- update_JxW_values | update_quadrature_points);
-
- const UpdateFlags local_face_flags (update_values);
-
- const UpdateFlags flags ( update_values | update_normal_vectors |
- update_quadrature_points |
- update_JxW_values);
-
- PerTaskData task_data (fe.dofs_per_cell,
- trace_reconstruct);
- ScratchData scratch (fe, fe_local,
- quadrature_formula,
- face_quadrature_formula,
- local_flags,
- local_face_flags,
- flags);
-
- WorkStream::run(dof_handler.begin_active(),
- dof_handler.end(),
- *this,
- &HDG<dim>::assemble_system_one_cell,
- &HDG<dim>::copy_local_to_global,
- scratch,
- task_data);
-}
+ const UpdateFlags local_face_flags (update_values);
+ const UpdateFlags flags ( update_values | update_normal_vectors |
+ update_quadrature_points |
+ update_JxW_values);
+ PerTaskData task_data (fe.dofs_per_cell,
+ trace_reconstruct);
+ ScratchData scratch (fe, fe_local,
+ quadrature_formula,
+ face_quadrature_formula,
+ local_flags,
+ local_face_flags,
+ flags);
+
+ WorkStream::run(dof_handler.begin_active(),
+ dof_handler.end(),
+ *this,
+ &HDG<dim>::assemble_system_one_cell,
+ &HDG<dim>::copy_local_to_global,
+ scratch,
+ task_data);
+ }
-// @sect4{HDG::assemble_system_one_cell}
-// The real work of the HDG program is done by @p assemble_system_one_cell.
-// Assembling the local matrices $A, B, C$ is done here, along with the
-// local contributions of the global matrix $D$.
-template <int dim>
-void
-HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
- ScratchData &scratch,
- PerTaskData &task_data)
-{
-// Construct iterator for dof_handler_local for FEValues reinit function.
- typename DoFHandler<dim>::active_cell_iterator
- loc_cell (&triangulation,
- cell->level(),
- cell->index(),
- &dof_handler_local);
- const unsigned int n_q_points = scratch.fe_values_local.get_quadrature().size();
- const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size();
- const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell;
+ // @sect4{HDG::assemble_system_one_cell}
+ // The real work of the HDG program is done by @p assemble_system_one_cell.
+ // Assembling the local matrices $A, B, C$ is done here, along with the
+ // local contributions of the global matrix $D$.
+ template <int dim>
+ void
+ HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ ScratchData &scratch,
+ PerTaskData &task_data)
+ {
+ // Construct iterator for dof_handler_local for FEValues reinit function.
+ typename DoFHandler<dim>::active_cell_iterator
+ loc_cell (&triangulation,
+ cell->level(),
+ cell->index(),
+ &dof_handler_local);
- const FEValuesExtractors::Vector fluxes (0);
- const FEValuesExtractors::Scalar scalar (dim);
+ const unsigned int n_q_points = scratch.fe_values_local.get_quadrature().size();
+ const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size();
- scratch.ll_matrix = 0;
- scratch.l_rhs = 0;
- if (!task_data.trace_reconstruct)
- {
- scratch.lf_matrix = 0;
- scratch.fl_matrix = 0;
- task_data.cell_matrix = 0;
- task_data.cell_vector = 0;
- }
- scratch.fe_values_local.reinit (loc_cell);
-
- // We first compute the cell-interior contribution to @p ll_matrix matrix
- // (referred to as matrix $A$ in the introduction) corresponding to
- // local-local coupling, as well as the local right-hand-side vector. We
- // store the values at each quadrature point for the basis functions, the
- // right-hand-side value, and the convection velocity.
- for (unsigned int q=0; q<n_q_points; ++q)
- {
- const double rhs_value
- = scratch.right_hand_side.value(scratch.fe_values_local.quadrature_point(q));
- const Tensor<1,dim> convection
- = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q));
- const double JxW = scratch.fe_values_local.JxW(q);
- for (unsigned int k=0; k<loc_dofs_per_cell; ++k)
- {
- scratch.q_phi[k] = scratch.fe_values_local[fluxes].value(k,q);
- scratch.q_phi_div[k] = scratch.fe_values_local[fluxes].divergence(k,q);
- scratch.u_phi[k] = scratch.fe_values_local[scalar].value(k,q);
- scratch.u_phi_grad[k] = scratch.fe_values_local[scalar].gradient(k,q);
- }
- for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<loc_dofs_per_cell; ++j)
- scratch.ll_matrix(i,j) += (
- scratch.q_phi[i] * scratch.q_phi[j]
- -
- scratch.q_phi_div[i] * scratch.u_phi[j]
- +
- scratch.u_phi[i] * scratch.q_phi_div[j]
- -
- (scratch.u_phi_grad[i] * convection) * scratch.u_phi[j]
- ) * JxW;
- scratch.l_rhs(i) += scratch.u_phi[i] * rhs_value * JxW;
- }
- }
+ const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell;
- // Face terms are assembled on all faces of all elements. This is in
- // contrast to more traditional DG methods, where each face is only visited
- // once in the assembly procedure.
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- {
- scratch.fe_face_values_local.reinit(loc_cell, face);
- scratch.fe_face_values.reinit(cell, face);
+ const FEValuesExtractors::Vector fluxes (0);
+ const FEValuesExtractors::Scalar scalar (dim);
- // The already obtained $\hat{u}$ values are needed when solving for the
- // local variables.
- if (task_data.trace_reconstruct)
- scratch.fe_face_values.get_function_values (solution, scratch.trace_values);
+ scratch.ll_matrix = 0;
+ scratch.l_rhs = 0;
+ if (!task_data.trace_reconstruct)
+ {
+ scratch.lf_matrix = 0;
+ scratch.fl_matrix = 0;
+ task_data.cell_matrix = 0;
+ task_data.cell_vector = 0;
+ }
+ scratch.fe_values_local.reinit (loc_cell);
+
+ // We first compute the cell-interior contribution to @p ll_matrix matrix
+ // (referred to as matrix $A$ in the introduction) corresponding to
+ // local-local coupling, as well as the local right-hand-side vector. We
+ // store the values at each quadrature point for the basis functions, the
+ // right-hand-side value, and the convection velocity, in order to have
+ // quick access to these fields.
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ const double rhs_value
+ = scratch.right_hand_side.value(scratch.fe_values_local.quadrature_point(q));
+ const Tensor<1,dim> convection
+ = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q));
+ const double JxW = scratch.fe_values_local.JxW(q);
+ for (unsigned int k=0; k<loc_dofs_per_cell; ++k)
+ {
+ scratch.q_phi[k] = scratch.fe_values_local[fluxes].value(k,q);
+ scratch.q_phi_div[k] = scratch.fe_values_local[fluxes].divergence(k,q);
+ scratch.u_phi[k] = scratch.fe_values_local[scalar].value(k,q);
+ scratch.u_phi_grad[k] = scratch.fe_values_local[scalar].gradient(k,q);
+ }
+ for (unsigned int i=0; i<loc_dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<loc_dofs_per_cell; ++j)
+ scratch.ll_matrix(i,j) += (
+ scratch.q_phi[i] * scratch.q_phi[j]
+ -
+ scratch.q_phi_div[i] * scratch.u_phi[j]
+ +
+ scratch.u_phi[i] * scratch.q_phi_div[j]
+ -
+ (scratch.u_phi_grad[i] * convection) * scratch.u_phi[j]
+ ) * JxW;
+ scratch.l_rhs(i) += scratch.u_phi[i] * rhs_value * JxW;
+ }
+ }
- for (unsigned int q=0; q<n_face_q_points; ++q)
- {
- const double JxW = scratch.fe_face_values.JxW(q);
- const Point<dim> quadrature_point =
- scratch.fe_face_values.quadrature_point(q);
- const Point<dim> normal = scratch.fe_face_values.normal_vector(q);
- const Tensor<1,dim> convection
- = scratch.convection_velocity.value(quadrature_point);
-
- // Here we compute the stabilization parameter discussed in the
- // introduction: since the diffusion is one and the diffusion length
- // scale is set to 1/5, it simply results in a contribution of 5 for
- // the diffusion part and the magnitude of convection through the
- // element boundary in a centered-like scheme for the convection
- // part.
- const double tau_stab = (5. +
- std::abs(convection * normal));
-
- // We store the non-zero flux and scalar values, making use of the
- // support_on_face information we calculated in @p ScratchData.
- for (unsigned int k=0; k<scratch.fe_local_support_on_face[face].size(); ++k)
- {
- const unsigned int kk=scratch.fe_local_support_on_face[face][k];
- scratch.q_phi[k] = scratch.fe_face_values_local[fluxes].value(kk,q);
- scratch.u_phi[k] = scratch.fe_face_values_local[scalar].value(kk,q);
- }
-
- // When @p trace_reconstruct=false, we are preparing assemble the
- // system for the skeleton variable $\lambda$. If this is the case,
- // we must assemble all local matrices associated with the problem:
- // local-local, local-face, face-local, and face-face. The
- // face-face matrix is stored as @p TaskData::cell_matrix, so that
- // it can be assembled into the global system by @p
- // copy_local_to_global.
- if (!task_data.trace_reconstruct)
- {
- for (unsigned int k=0; k<scratch.fe_support_on_face[face].size(); ++k)
- scratch.tr_phi[k] =
- scratch.fe_face_values.shape_value(scratch.fe_support_on_face[face][k],q);
- for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
- for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
- {
- const unsigned int ii=scratch.fe_local_support_on_face[face][i];
- const unsigned int jj=scratch.fe_support_on_face[face][j];
- scratch.lf_matrix(ii,jj) += (
- (scratch.q_phi[i] * normal
- +
- (convection * normal -
- tau_stab) * scratch.u_phi[i])
- * scratch.tr_phi[j]
- ) * JxW;
-
- // Note the sign of the face-local matrix. We negate the
- // sign during assembly here so that we can use the
- // FullMatrix::mmult with addition when computing the
- // Schur complement.
- scratch.fl_matrix(jj,ii) -= (
- (scratch.q_phi[i] * normal
- +
- tau_stab * scratch.u_phi[i])
- * scratch.tr_phi[j]
- ) * JxW;
- }
+ // Face terms are assembled on all faces of all elements. This is in
+ // contrast to more traditional DG methods, where each face is only visited
+ // once in the assembly procedure.
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ {
+ scratch.fe_face_values_local.reinit(loc_cell, face);
+ scratch.fe_face_values.reinit(cell, face);
+
+ // The already obtained $\hat{u}$ values are needed when solving for the
+ // local variables.
+ if (task_data.trace_reconstruct)
+ scratch.fe_face_values.get_function_values (solution, scratch.trace_values);
+
+ for (unsigned int q=0; q<n_face_q_points; ++q)
+ {
+ const double JxW = scratch.fe_face_values.JxW(q);
+ const Point<dim> quadrature_point =
+ scratch.fe_face_values.quadrature_point(q);
+ const Point<dim> normal = scratch.fe_face_values.normal_vector(q);
+ const Tensor<1,dim> convection
+ = scratch.convection_velocity.value(quadrature_point);
+
+ // Here we compute the stabilization parameter discussed in the
+ // introduction: since the diffusion is one and the diffusion
+ // length scale is set to 1/5, it simply results in a contribution
+ // of 5 for the diffusion part and the magnitude of convection
+ // through the element boundary in a centered scheme for the
+ // convection part.
+ const double tau_stab = (5. +
+ std::abs(convection * normal));
+
+ // We store the non-zero flux and scalar values, making use of the
+ // support_on_face information we created in @p ScratchData.
+ for (unsigned int k=0; k<scratch.fe_local_support_on_face[face].size(); ++k)
+ {
+ const unsigned int kk=scratch.fe_local_support_on_face[face][k];
+ scratch.q_phi[k] = scratch.fe_face_values_local[fluxes].value(kk,q);
+ scratch.u_phi[k] = scratch.fe_face_values_local[scalar].value(kk,q);
+ }
- for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
- for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
- {
- const unsigned int ii=scratch.fe_support_on_face[face][i];
- const unsigned int jj=scratch.fe_support_on_face[face][j];
- task_data.cell_matrix(ii,jj) += (
- (convection * normal - tau_stab) *
- scratch.tr_phi[i] * scratch.tr_phi[j]
- ) * JxW;
- }
+ // When @p trace_reconstruct=false, we are preparing to assemble the
+ // system for the skeleton variable $\hat{u}$. If this is the case,
+ // we must assemble all local matrices associated with the problem:
+ // local-local, local-face, face-local, and face-face. The
+ // face-face matrix is stored as @p TaskData::cell_matrix, so that
+ // it can be assembled into the global system by @p
+ // copy_local_to_global.
+ if (!task_data.trace_reconstruct)
+ {
+ for (unsigned int k=0; k<scratch.fe_support_on_face[face].size(); ++k)
+ scratch.tr_phi[k] =
+ scratch.fe_face_values.shape_value(scratch.fe_support_on_face[face][k],q);
+ for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+ for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ const unsigned int jj=scratch.fe_support_on_face[face][j];
+ scratch.lf_matrix(ii,jj) += (
+ (scratch.q_phi[i] * normal
+ +
+ (convection * normal -
+ tau_stab) * scratch.u_phi[i])
+ * scratch.tr_phi[j]
+ ) * JxW;
+
+ // Note the sign of the face-local matrix. We negate the
+ // sign during assembly here so that we can use the
+ // FullMatrix::mmult with addition when computing the
+ // Schur complement.
+ scratch.fl_matrix(jj,ii) -= (
+ (scratch.q_phi[i] * normal
+ +
+ tau_stab * scratch.u_phi[i])
+ * scratch.tr_phi[j]
+ ) * JxW;
+ }
- if (cell->face(face)->at_boundary()
- &&
- (cell->face(face)->boundary_indicator() == 1))
- {
- const double neumann_value =
- - scratch.exact_solution.gradient (quadrature_point) * normal
- + convection * normal * scratch.exact_solution.value(quadrature_point);
- for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+ for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+ for (unsigned int j=0; j<scratch.fe_support_on_face[face].size(); ++j)
{
const unsigned int ii=scratch.fe_support_on_face[face][i];
- task_data.cell_vector(ii) += scratch.tr_phi[i] * neumann_value * JxW;
+ const unsigned int jj=scratch.fe_support_on_face[face][j];
+ task_data.cell_matrix(ii,jj) += (
+ (convection * normal - tau_stab) *
+ scratch.tr_phi[i] * scratch.tr_phi[j]
+ ) * JxW;
}
- }
- }
- // This last term adds the contribution of the term $\left<w,\tau
- // u_h\right>_{\partial \mathcal T}$ to the local matrix. As opposed
- // to the face matrices above, we need it in both assembly stages.
- for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
- for (unsigned int j=0; j<scratch.fe_local_support_on_face[face].size(); ++j)
- {
- const unsigned int ii=scratch.fe_local_support_on_face[face][i];
- const unsigned int jj=scratch.fe_local_support_on_face[face][j];
- scratch.ll_matrix(ii,jj) += tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW;
+ if (cell->face(face)->at_boundary()
+ &&
+ (cell->face(face)->boundary_indicator() == 1))
+ {
+ const double neumann_value =
+ - scratch.exact_solution.gradient (quadrature_point) * normal
+ + convection * normal * scratch.exact_solution.value(quadrature_point);
+ for (unsigned int i=0; i<scratch.fe_support_on_face[face].size(); ++i)
+ {
+ const unsigned int ii=scratch.fe_support_on_face[face][i];
+ task_data.cell_vector(ii) += scratch.tr_phi[i] * neumann_value * JxW;
+ }
+ }
}
- // When @p trace_reconstruct=true, we are solving for the local
- // solutions on an element by element basis. The local
- // right-hand-side is calculated by replacing the basis functions @p
- // tr_phi in the @p lf_matrix computation by the computed values @p
- // trace_values. Of course, the sign of the matrix is now minus
- // since we have moved everything to the other side of the equation.
- if (task_data.trace_reconstruct)
+ // This last term adds the contribution of the term $\left<w,\tau
+ // u_h\right>_{\partial \mathcal T}$ to the local matrix. As opposed
+ // to the face matrices above, we need it in both assembly stages.
for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
- {
- const unsigned int ii=scratch.fe_local_support_on_face[face][i];
- scratch.l_rhs(ii) -= (scratch.q_phi[i] * normal
- +
- scratch.u_phi[i] * (convection * normal - tau_stab)
- ) * scratch.trace_values[q] * JxW;
- }
- }
- }
+ for (unsigned int j=0; j<scratch.fe_local_support_on_face[face].size(); ++j)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ const unsigned int jj=scratch.fe_local_support_on_face[face][j];
+ scratch.ll_matrix(ii,jj) += tau_stab * scratch.u_phi[i] * scratch.u_phi[j] * JxW;
+ }
+
+ // When @p trace_reconstruct=true, we are solving for the local
+ // solutions on an element by element basis. The local
+ // right-hand-side is calculated by replacing the basis functions @p
+ // tr_phi in the @p lf_matrix computation by the computed values @p
+ // trace_values. Of course, the sign of the matrix is now minus
+ // since we have moved everything to the other side of the equation.
+ if (task_data.trace_reconstruct)
+ for (unsigned int i=0; i<scratch.fe_local_support_on_face[face].size(); ++i)
+ {
+ const unsigned int ii=scratch.fe_local_support_on_face[face][i];
+ scratch.l_rhs(ii) -= (scratch.q_phi[i] * normal
+ +
+ scratch.u_phi[i] * (convection * normal - tau_stab)
+ ) * scratch.trace_values[q] * JxW;
+ }
+ }
+ }
- // Once assembly of all of the local contributions is complete, we must either:
- // (1) assemble the global system, or (2) compute the local solution values and
- // save them.
- // In either case, the first step is to invert the local-local matrix.
- scratch.ll_matrix.gauss_jordan();
+ // Once assembly of all of the local contributions is complete, we must either:
+ // (1) assemble the global system, or (2) compute the local solution values and
+ // save them.
+ // In either case, the first step is to invert the local-local matrix.
+ scratch.ll_matrix.gauss_jordan();
- // For (1), we compute the Schur complement and add it to the @p
- // cell_matrix, matrix $D$ in the introduction.
- if (task_data.trace_reconstruct == false)
- {
- scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix);
- scratch.tmp_matrix.vmult_add(task_data.cell_vector, scratch.l_rhs);
- scratch.tmp_matrix.mmult(task_data.cell_matrix, scratch.lf_matrix, true);
- cell->get_dof_indices(task_data.dof_indices);
- }
- // For (2), we are simply solving (ll_matrix).(solution_local) = (l_rhs).
- // Hence, we multiply @p l_rhs by our already inverted local-local matrix
- // and store the result using the <code>set_dof_values</code> function.
- else
- {
- scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs);
- loc_cell->set_dof_values(scratch.tmp_rhs, solution_local);
- }
-}
+ // For (1), we compute the Schur complement and add it to the @p
+ // cell_matrix, matrix $D$ in the introduction.
+ if (task_data.trace_reconstruct == false)
+ {
+ scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix);
+ scratch.tmp_matrix.vmult_add(task_data.cell_vector, scratch.l_rhs);
+ scratch.tmp_matrix.mmult(task_data.cell_matrix, scratch.lf_matrix, true);
+ cell->get_dof_indices(task_data.dof_indices);
+ }
+ // For (2), we are simply solving (ll_matrix).(solution_local) = (l_rhs).
+ // Hence, we multiply @p l_rhs by our already inverted local-local matrix
+ // and store the result using the <code>set_dof_values</code> function.
+ else
+ {
+ scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs);
+ loc_cell->set_dof_values(scratch.tmp_rhs, solution_local);
+ }
+ }
-// @sect4{HDG::solve}
-// The skeleton solution is solved for by using a BiCGStab solver with
-// identity preconditioner.
-template <int dim>
-void HDG<dim>::solve ()
-{
- SolverControl solver_control (system_matrix.m()*10,
- 1e-11*system_rhs.l2_norm());
- SolverBicgstab<> solver (solver_control, false);
- solver.solve (system_matrix, solution, system_rhs,
- PreconditionIdentity());
+ // @sect4{HDG::copy_local_to_global}
+ // If we are in the first step of the solution, i.e. @p trace_reconstruct=false,
+ // then we assemble the local matrices into the global system.
+ template <int dim>
+ void HDG<dim>::copy_local_to_global(const PerTaskData &data)
+ {
+ if (data.trace_reconstruct == false)
+ constraints.distribute_local_to_global (data.cell_matrix,
+ data.cell_vector,
+ data.dof_indices,
+ system_matrix, system_rhs);
+ }
- std::cout << " Number of BiCGStab iterations: " << solver_control.last_step()
- << std::endl;
- system_matrix.clear();
- sparsity_pattern.reinit(0,0,0,1);
- constraints.distribute(solution);
+ // @sect4{HDG::solve}
+ // The skeleton solution is solved for by using a BiCGStab solver with
+ // identity preconditioner.
+ template <int dim>
+ void HDG<dim>::solve ()
+ {
+ SolverControl solver_control (system_matrix.m()*10,
+ 1e-11*system_rhs.l2_norm());
+ SolverBicgstab<> solver (solver_control, false);
+ solver.solve (system_matrix, solution, system_rhs,
+ PreconditionIdentity());
- // Once we have solved for the skeleton solution,
- // we can solve for the local solutions in an element-by-element
- // fashion. We do this by re-using the same @p assemble_system function
- // but switching @p trace_reconstruct to true.
- assemble_system(true);
-}
+ std::cout << " Number of BiCGStab iterations: " << solver_control.last_step()
+ << std::endl;
+
+ system_matrix.clear();
+ sparsity_pattern.reinit(0,0,0,1);
+
+ constraints.distribute(solution);
+
+ // Once we have solved for the skeleton solution,
+ // we can solve for the local solutions in an element-by-element
+ // fashion. We do this by re-using the same @p assemble_system function
+ // but switching @p trace_reconstruct to true.
+ assemble_system(true);
+ }
// SolutionAndGradient class introduced above that contains the analytic
// parts of either of them. Eventually, we also compute the L2-error of the
// post-processed solution and add the results into the convergence table.
-template <int dim>
-void
-HDG<dim>::postprocess()
-{
+ template <int dim>
+ void
+ HDG<dim>::postprocess()
{
- const QGauss<dim> quadrature_formula(fe_u_post.degree+1);
- const UpdateFlags local_flags (update_values);
- const UpdateFlags flags ( update_values | update_gradients |
- update_JxW_values);
+ {
+ const QGauss<dim> quadrature_formula(fe_u_post.degree+1);
+ const UpdateFlags local_flags (update_values);
+ const UpdateFlags flags ( update_values | update_gradients |
+ update_JxW_values);
+
+ PostProcessScratchData scratch (fe_u_post, fe_local,
+ quadrature_formula,
+ local_flags,
+ flags);
+
+ WorkStream::run(dof_handler_u_post.begin_active(),
+ dof_handler_u_post.end(),
+ std_cxx1x::bind (&HDG<dim>::postprocess_one_cell,
+ std_cxx1x::ref(*this),
+ std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
+ std_cxx1x::function<void(const unsigned int &)>(),
+ scratch,
+ 0U);
+ }
- PostProcessScratchData scratch (fe_u_post, fe_local,
- quadrature_formula,
- local_flags,
- flags);
-
- WorkStream::run(dof_handler_u_post.begin_active(),
- dof_handler_u_post.end(),
- std_cxx1x::bind (&HDG<dim>::postprocess_one_cell,
- std_cxx1x::ref(*this),
- std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
- std_cxx1x::function<void(const unsigned int&)>(),
- scratch,
- 0U);
+ Vector<float> difference_per_cell (triangulation.n_active_cells());
+
+ ComponentSelectFunction<dim> value_select (dim, dim+1);
+ VectorTools::integrate_difference (dof_handler_local,
+ solution_local,
+ SolutionAndGradient<dim>(),
+ difference_per_cell,
+ QGauss<dim>(fe.degree+2),
+ VectorTools::L2_norm,
+ &value_select);
+ const double L2_error = difference_per_cell.l2_norm();
+
+ ComponentSelectFunction<dim> gradient_select (std::pair<unsigned int,unsigned int>(0, dim),
+ dim+1);
+ VectorTools::integrate_difference (dof_handler_local,
+ solution_local,
+ SolutionAndGradient<dim>(),
+ difference_per_cell,
+ QGauss<dim>(fe.degree+2),
+ VectorTools::L2_norm,
+ &gradient_select);
+ const double grad_error = difference_per_cell.l2_norm();
+
+ VectorTools::integrate_difference (dof_handler_u_post,
+ solution_u_post,
+ Solution<dim>(),
+ difference_per_cell,
+ QGauss<dim>(fe.degree+3),
+ VectorTools::L2_norm);
+ const double post_error = difference_per_cell.l2_norm();
+
+ convergence_table.add_value("cells", triangulation.n_active_cells());
+ convergence_table.add_value("dofs", dof_handler.n_dofs());
+ convergence_table.add_value("val L2", L2_error);
+ convergence_table.add_value("grad L2", grad_error);
+ convergence_table.add_value("val L2-post", post_error);
}
- Vector<float> difference_per_cell (triangulation.n_active_cells());
-
- ComponentSelectFunction<dim> value_select (dim, dim+1);
- VectorTools::integrate_difference (dof_handler_local,
- solution_local,
- SolutionAndGradient<dim>(),
- difference_per_cell,
- QGauss<dim>(fe.degree+2),
- VectorTools::L2_norm,
- &value_select);
- const double L2_error = difference_per_cell.l2_norm();
-
- ComponentSelectFunction<dim> gradient_select (std::pair<unsigned int,unsigned int>(0, dim),
- dim+1);
- VectorTools::integrate_difference (dof_handler_local,
- solution_local,
- SolutionAndGradient<dim>(),
- difference_per_cell,
- QGauss<dim>(fe.degree+2),
- VectorTools::L2_norm,
- &gradient_select);
- const double grad_error = difference_per_cell.l2_norm();
-
- VectorTools::integrate_difference (dof_handler_u_post,
- solution_u_post,
- Solution<dim>(),
- difference_per_cell,
- QGauss<dim>(fe.degree+3),
- VectorTools::L2_norm);
- const double post_error = difference_per_cell.l2_norm();
-
- convergence_table.add_value("cells", triangulation.n_active_cells());
- convergence_table.add_value("dofs", dof_handler.n_dofs());
- convergence_table.add_value("val L2", L2_error);
- convergence_table.add_value("grad L2", grad_error);
- convergence_table.add_value("val L2-post", post_error);
-}
-
// @sect4{HDG::postprocess_one_cell}
// discussion in the introduction, we need to set up a system that projects
// the gradient part of the DG solution onto the gradient of the
// post-processed variable. Moreover, we need to set the average of the new
- // post-processed variable to be equal the average of the scalar DG solution
+ // post-processed variable to equal the average of the scalar DG solution
// on the cell.
//
// More technically speaking, the projection of the gradient is a system
// row would correspond to the constant part already and deleting e.g. the
// last row would give us a singular system. This way, our program can also
// be used for those elements.
-template <int dim>
-void
-HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
- PostProcessScratchData &scratch,
- unsigned int &)
-{
- typename DoFHandler<dim>::active_cell_iterator
- loc_cell (&triangulation,
- cell->level(),
- cell->index(),
- &dof_handler_local);
+ template <int dim>
+ void
+ HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+ PostProcessScratchData &scratch,
+ unsigned int &)
+ {
+ typename DoFHandler<dim>::active_cell_iterator
+ loc_cell (&triangulation,
+ cell->level(),
+ cell->index(),
+ &dof_handler_local);
- scratch.fe_values_local.reinit (loc_cell);
- scratch.fe_values.reinit(cell);
+ scratch.fe_values_local.reinit (loc_cell);
+ scratch.fe_values.reinit(cell);
- FEValuesExtractors::Vector fluxes(0);
- FEValuesExtractors::Scalar scalar(dim);
+ FEValuesExtractors::Vector fluxes(0);
+ FEValuesExtractors::Scalar scalar(dim);
- const unsigned int n_q_points = scratch.fe_values.get_quadrature().size();
- const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell;
+ const unsigned int n_q_points = scratch.fe_values.get_quadrature().size();
+ const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell;
- scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values);
- scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients);
+ scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values);
+ scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients);
- double sum = 0;
- for (unsigned int i=1; i<dofs_per_cell; ++i)
- {
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- {
- sum = 0;
- for (unsigned int q=0; q<n_q_points; ++q)
- sum += (scratch.fe_values.shape_grad(i,q) *
- scratch.fe_values.shape_grad(j,q)
- ) * scratch.fe_values.JxW(q);
- scratch.cell_matrix(i,j) = sum;
- }
-
- sum = 0;
- for (unsigned int q=0; q<n_q_points; ++q)
- sum -= (scratch.fe_values.shape_grad(i,q) * scratch.u_gradients[q]
- ) * scratch.fe_values.JxW(q);
- scratch.cell_rhs(i) = sum;
- }
- for (unsigned int j=0; j<dofs_per_cell; ++j)
+ double sum = 0;
+ for (unsigned int i=1; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum += (scratch.fe_values.shape_grad(i,q) *
+ scratch.fe_values.shape_grad(j,q)
+ ) * scratch.fe_values.JxW(q);
+ scratch.cell_matrix(i,j) = sum;
+ }
+
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum -= (scratch.fe_values.shape_grad(i,q) * scratch.u_gradients[q]
+ ) * scratch.fe_values.JxW(q);
+ scratch.cell_rhs(i) = sum;
+ }
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ sum = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ sum += scratch.fe_values.shape_value(j,q) * scratch.fe_values.JxW(q);
+ scratch.cell_matrix(0,j) = sum;
+ }
{
sum = 0;
for (unsigned int q=0; q<n_q_points; ++q)
- sum += scratch.fe_values.shape_value(j,q) * scratch.fe_values.JxW(q);
- scratch.cell_matrix(0,j) = sum;
+ sum += scratch.u_values[q] * scratch.fe_values.JxW(q);
+ scratch.cell_rhs(0) = sum;
}
- {
- sum = 0;
- for (unsigned int q=0; q<n_q_points; ++q)
- sum += scratch.u_values[q] * scratch.fe_values.JxW(q);
- scratch.cell_rhs(0) = sum;
- }
- // Having assembled all terms, we can again go on and solve the linear
- // system. We again invert the matrix and then multiply the inverse by the
- // right hand side. An alternative (and more numerically stable) would have
- // been to only factorize the matrix and apply the factorization.
- scratch.cell_matrix.gauss_jordan();
- scratch.cell_matrix.vmult(scratch.cell_sol, scratch.cell_rhs);
- cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
-}
+ // Having assembled all terms, we can again go on and solve the linear
+ // system. We invert the matrix and then multiply the inverse by the
+ // right hand side. An alternative (and more numerically stable) method would have
+ // been to only factorize the matrix and apply the factorization.
+ scratch.cell_matrix.gauss_jordan();
+ scratch.cell_matrix.vmult(scratch.cell_sol, scratch.cell_rhs);
+ cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
+ }
-// @sect4{HDG::output_results}
-// We have 3 sets of results that we would like to output: the local solution,
-// the post-processed local solution, and the skeleton solution. The former 2
-// both `live' on element volumes, wheras the latter lives on codimention-1 surfaces
-// of the triangulation. Our @p output_results function writes all local solutions
-// to the same vtk file, even though they correspond to different <code>DoFHandler</code>
-// objects. The graphical output for the skeleton variable is done through
-// use of the <code>DataOutFaces</code> class.
-template <int dim>
-void HDG<dim>::output_results (const unsigned int cycle)
-{
- std::string filename;
- switch (refinement_mode)
- {
- case global_refinement:
- filename = "solution-global";
- break;
- case adaptive_refinement:
- filename = "solution-adaptive";
- break;
- default:
- Assert (false, ExcNotImplemented());
- }
+ // @sect4{HDG::output_results}
+ // We have 3 sets of results that we would like to output: the local solution,
+ // the post-processed local solution, and the skeleton solution. The former 2
+ // both 'live' on element volumes, wheras the latter lives on codimention-1 surfaces
+ // of the triangulation. Our @p output_results function writes all local solutions
+ // to the same vtk file, even though they correspond to different <code>DoFHandler</code>
+ // objects. The graphical output for the skeleton variable is done through
+ // use of the <code>DataOutFaces</code> class.
+ template <int dim>
+ void HDG<dim>::output_results (const unsigned int cycle)
+ {
+ std::string filename;
+ switch (refinement_mode)
+ {
+ case global_refinement:
+ filename = "solution-global";
+ break;
+ case adaptive_refinement:
+ filename = "solution-adaptive";
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
- std::string face_out(filename);
- face_out += "-face";
-
- filename += "-q" + Utilities::int_to_string(fe.degree,1);
- filename += "-" + Utilities::int_to_string(cycle,2);
- filename += ".vtk";
- std::ofstream output (filename.c_str());
-
- DataOut<dim> data_out;
-
-// We first define the names and types of the local solution,
-// and add the data to @p data_out.
- std::vector<std::string> names (dim, "gradient");
- names.push_back ("solution");
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- component_interpretation
- (dim+1, DataComponentInterpretation::component_is_part_of_vector);
- component_interpretation[dim]
- = DataComponentInterpretation::component_is_scalar;
- data_out.add_data_vector (dof_handler_local, solution_local,
- names, component_interpretation);
-
-// The second data item we add is the post-processed solution.
-// In this case, it is a single scalar variable belonging to
-// a different DoFHandler.
- std::vector<std::string> post_name(1,"u_post");
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- post_comp_type(1, DataComponentInterpretation::component_is_scalar);
- data_out.add_data_vector (dof_handler_u_post, solution_u_post,
- post_name, post_comp_type);
-
- data_out.build_patches (fe.degree);
- data_out.write_vtk (output);
-
- face_out += "-q" + Utilities::int_to_string(fe.degree,1);
- face_out += "-" + Utilities::int_to_string(cycle,2);
- face_out += ".vtk";
- std::ofstream face_output (face_out.c_str());
+ std::string face_out(filename);
+ face_out += "-face";
+
+ filename += "-q" + Utilities::int_to_string(fe.degree,1);
+ filename += "-" + Utilities::int_to_string(cycle,2);
+ filename += ".vtk";
+ std::ofstream output (filename.c_str());
+
+ DataOut<dim> data_out;
+
+ // We first define the names and types of the local solution,
+ // and add the data to @p data_out.
+ std::vector<std::string> names (dim, "gradient");
+ names.push_back ("solution");
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ component_interpretation
+ (dim+1, DataComponentInterpretation::component_is_part_of_vector);
+ component_interpretation[dim]
+ = DataComponentInterpretation::component_is_scalar;
+ data_out.add_data_vector (dof_handler_local, solution_local,
+ names, component_interpretation);
+
+ // The second data item we add is the post-processed solution.
+ // In this case, it is a single scalar variable belonging to
+ // a different DoFHandler.
+ std::vector<std::string> post_name(1,"u_post");
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ post_comp_type(1, DataComponentInterpretation::component_is_scalar);
+ data_out.add_data_vector (dof_handler_u_post, solution_u_post,
+ post_name, post_comp_type);
+
+ data_out.build_patches (fe.degree);
+ data_out.write_vtk (output);
+
+ face_out += "-q" + Utilities::int_to_string(fe.degree,1);
+ face_out += "-" + Utilities::int_to_string(cycle,2);
+ face_out += ".vtk";
+ std::ofstream face_output (face_out.c_str());
// The <code>DataOutFaces</code> class works analagously to the <code>DataOut</code>
// class when we have a <code>DoFHandler</code> that defines the solution on
// the skeleton of the triangulation. We treat it as such here, and the code is
// similar to that above.
- DataOutFaces<dim> data_out_face(false);
- std::vector<std::string> face_name(1,"lambda");
- std::vector<DataComponentInterpretation::DataComponentInterpretation>
- face_component_type(1, DataComponentInterpretation::component_is_scalar);
-
- data_out_face.add_data_vector (dof_handler,
- solution,
- face_name,
- face_component_type);
-
- data_out_face.build_patches (fe.degree);
- data_out_face.write_vtk (face_output);
-}
+ DataOutFaces<dim> data_out_face(false);
+ std::vector<std::string> face_name(1,"u_hat");
+ std::vector<DataComponentInterpretation::DataComponentInterpretation>
+ face_component_type(1, DataComponentInterpretation::component_is_scalar);
+
+ data_out_face.add_data_vector (dof_handler,
+ solution,
+ face_name,
+ face_component_type);
+
+ data_out_face.build_patches (fe.degree);
+ data_out_face.write_vtk (face_output);
+ }
// @sect4{HDG::refine_grid}
// The adaptive_refinement mode uses the <code>KellyErrorEstimator</code> to
// give a decent indication of the non-regular regions in the scalar local
// solutions.
-template <int dim>
-void HDG<dim>::refine_grid (const unsigned int cycle)
-{
- if (cycle == 0)
- {
- GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1);
- triangulation.refine_global(3-dim);
- }
- else
- switch (refinement_mode)
- {
- case global_refinement:
+ template <int dim>
+ void HDG<dim>::refine_grid (const unsigned int cycle)
+ {
+ if (cycle == 0)
{
- triangulation.clear();
- GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
- triangulation.refine_global(3-dim+cycle/2);
- break;
+ GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1);
+ triangulation.refine_global(3-dim);
}
+ else
+ switch (refinement_mode)
+ {
+ case global_refinement:
+ {
+ triangulation.clear();
+ GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1);
+ triangulation.refine_global(3-dim+cycle/2);
+ break;
+ }
- case adaptive_refinement:
- {
- Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+ case adaptive_refinement:
+ {
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
- FEValuesExtractors::Scalar scalar(dim);
- typename FunctionMap<dim>::type neumann_boundary;
- KellyErrorEstimator<dim>::estimate (dof_handler_local,
- QGauss<dim-1>(3),
- neumann_boundary,
- solution_local,
- estimated_error_per_cell,
- fe_local.component_mask(scalar));
+ FEValuesExtractors::Scalar scalar(dim);
+ typename FunctionMap<dim>::type neumann_boundary;
+ KellyErrorEstimator<dim>::estimate (dof_handler_local,
+ QGauss<dim-1>(3),
+ neumann_boundary,
+ solution_local,
+ estimated_error_per_cell,
+ fe_local.component_mask(scalar));
- GridRefinement::refine_and_coarsen_fixed_number (triangulation,
- estimated_error_per_cell,
- 0.3, 0.);
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.);
- triangulation.execute_coarsening_and_refinement ();
+ triangulation.execute_coarsening_and_refinement ();
- break;
- }
+ break;
+ }
- default:
+ default:
+ {
+ Assert (false, ExcNotImplemented());
+ }
+ }
+
+ // Just as in step-7, we set the boundary indicator of two of the faces to 1
+ // where we want to specify Neumann boundary conditions instead of Dirichlet
+ // conditions. Since we re-create the triangulation every time for global
+ // refinement, the flags are set in every refinement step, not just at the
+ // beginning.
+ typename Triangulation<dim>::cell_iterator
+ cell = triangulation.begin (),
+ endc = triangulation.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->at_boundary())
+ if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12)
+ ||
+ (std::fabs(cell->face(face)->center()(1) - (-1)) < 1e-12))
+ cell->face(face)->set_boundary_indicator (1);
+ }
+
+ // @sect4{HDG::run}
+ // The functionality here is basically the same as <code>Step-7</code>.
+ // We loop over 10 cycles, refining the grid on each one. At the end,
+ // convergence tables are created.
+ template <int dim>
+ void HDG<dim>::run ()
+ {
+ for (unsigned int cycle=0; cycle<10; ++cycle)
{
- Assert (false, ExcNotImplemented());
- }
+ std::cout << "Cycle " << cycle << ':' << std::endl;
+
+ refine_grid (cycle);
+ setup_system ();
+ assemble_system (false);
+ solve ();
+ postprocess();
+ output_results (cycle);
}
- // Just as in step-7, we set the boundary indicator of one of the faces to 1
- // where we want to specify Neumann boundary conditions instead of Dirichlet
- // conditions. Since we re-create the triangulation every time for global
- // refinement, the flags are set in every refinement step, not just at the
- // beginning.
- typename Triangulation<dim>::cell_iterator
- cell = triangulation.begin (),
- endc = triangulation.end();
- for (; cell!=endc; ++cell)
- for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- if (cell->face(face)->at_boundary())
- if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12)
- ||
- (std::fabs(cell->face(face)->center()(1) - (-1)) < 1e-12))
- cell->face(face)->set_boundary_indicator (1);
-}
-// @sect4{HDG::run}
-// The functionality here is basically the same as <code>Step-7</code>.
-// We loop over 10 cycles, refining the grid on each one. At the end,
-// convergence tables are created.
-template <int dim>
-void HDG<dim>::run ()
-{
- for (unsigned int cycle=0; cycle<10; ++cycle)
- {
- std::cout << "Cycle " << cycle << ':' << std::endl;
-
- refine_grid (cycle);
- setup_system ();
- assemble_system (false);
- solve ();
- postprocess();
- output_results (cycle);
- }
+ convergence_table.set_precision("val L2", 3);
+ convergence_table.set_scientific("val L2", true);
+ convergence_table.set_precision("grad L2", 3);
+ convergence_table.set_scientific("grad L2", true);
+ convergence_table.set_precision("val L2-post", 3);
+ convergence_table.set_scientific("val L2-post", true);
+ // There is one minor change for the convergence table compared to step-7:
+ // Since we did not refine our mesh by a factor two in each cycle (but
+ // rather used the sequence 2, 3, 4, 6, 8, 12, ...), we need to tell the
+ // convergence rate evaluation about this. We do this by setting the
+ // number of cells as a reference column and additionally specifying the
+ // dimension of the problem, which gives the necessary information for the
+ // relation between number of cells and mesh size.
+ if (refinement_mode == global_refinement)
+ {
+ convergence_table
+ .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
+ convergence_table
+ .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
+ convergence_table
+ .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim);
+ }
+ convergence_table.write_text(std::cout);
+ }
- convergence_table.set_precision("val L2", 3);
- convergence_table.set_scientific("val L2", true);
- convergence_table.set_precision("grad L2", 3);
- convergence_table.set_scientific("grad L2", true);
- convergence_table.set_precision("val L2-post", 3);
- convergence_table.set_scientific("val L2-post", true);
+} // end of namespace Step51
- // There is one minor change for the convergence table compared to step-7:
- // Since we did not refine our mesh by a factor two in each cycle (but
- // rather used the sequence 2, 3, 4, 6, 8, 12, ...), we need to tell the
- // convergence rate evaluation about this. We do this by setting the number
- // of cells as a reference column and additionally specifying the dimension
- // of the problem, which gives the computation the necessary information for
- // how much the mesh was refinement given a certain increase in the number
- // of cells.
- if (refinement_mode == global_refinement)
- {
- convergence_table
- .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
- convergence_table
- .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim);
- convergence_table
- .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim);
- }
- convergence_table.write_text(std::cout);
-}
-} // end of namespace Step51
int main (int argc, char **argv)
{
--- /dev/null
+
+// ---------------------------------------------------------------------
+// $Id: graph_coloring.h 30494 2013-08-26 10:04:44Z kronbichler $
+//
+// Copyright (C) 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__graph_coloring_h
+#define __deal2__graph_coloring_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <boost/unordered_map.hpp>
+#include <boost/unordered_set.hpp>
+
+#include <vector>
+
+DEAL_II_NAMESPACE_OPEN
+
+/// This namespace contains the functions necessary to color a graph.
+namespace graph_coloring {
+
+ /**
+ * Create the partitioning using a simplified version of the Cuthill-McKee
+ * algorithm (Breadth First Search algorithm).
+ */
+ template <typename Iterator>
+ std::vector<std::vector<Iterator> > create_partitioning(Iterator const &begin,
+ typename identity<Iterator>::type const &end,
+ std::vector<types::global_dof_index> (*get_conflict_indices) (Iterator &))
+ {
+ std::vector<std::vector<Iterator> > partitioning(1,std::vector<Iterator> (1,begin));
+ // Number of iterators.
+ unsigned int n_iterators(0);
+ // Create a map from conflict indices to iterators
+ boost::unordered_map<types::global_dof_index,std::vector<Iterator> > indices_to_iterators;
+ for (Iterator it=begin; it!=end; ++it)
+ {
+ std::vector<types::global_dof_index> conflict_indices = (*get_conflict_indices)(it);
+ const unsigned int n_conflict_indices(conflict_indices.size());
+ for (unsigned int i=0; i<n_conflict_indices; ++i)
+ indices_to_iterators[conflict_indices[i]].push_back(it);
+ ++n_iterators;
+ }
+
+ // Create the partitioning.
+ std::set<Iterator> used_it;
+ used_it.insert(begin);
+ while (used_it.size()!=n_iterators)
+ {
+ typename std::vector<Iterator>::iterator vector_it(partitioning.back().begin());
+ typename std::vector<Iterator>::iterator vector_end(partitioning.back().end());
+ std::vector<Iterator> new_zone;
+ for (; vector_it!=vector_end; ++vector_it)
+ {
+ std::vector<types::global_dof_index> conflict_indices = (*get_conflict_indices)(*vector_it);
+ const unsigned int n_conflict_indices(conflict_indices.size());
+ for (unsigned int i=0; i<n_conflict_indices; ++i)
+ {
+ std::vector<Iterator> iterator_vector(indices_to_iterators[conflict_indices[i]]);
+ for (unsigned int j=0; j<iterator_vector.size(); ++j)
+ {
+ // Check that the iterator is not associated to a zone yet.
+ if (used_it.count(iterator_vector[j])==0)
+ {
+ new_zone.push_back(iterator_vector[j]);
+ used_it.insert(iterator_vector[j]);
+ }
+ }
+ }
+ }
+ // If there are iterators in the new zone, then the zone is added to the
+ // partition. Otherwise, the graph is disconnected and we need to find
+ // an iterator on the other part of the graph.
+ if (new_zone.size()!=0)
+ partitioning.push_back(new_zone);
+ else
+ for (Iterator it=begin; it!=end; ++it)
+ if (used_it.count(it)==0)
+ {
+ partitioning.push_back(std::vector<Iterator> (1,it));
+ break;
+ }
+ }
+
+ return partitioning;
+ }
+
+
+
+ /**
+ * This function uses DSATUR (Degree SATURation) to color one zone of the
+ * partition. DSATUR works as follows:
+ * -# Arrange the vertices by decreasing order of degrees.
+ * -# Color a vertex of maximal degree with color 1.
+ * -# Choose a vertex with a maximal saturation degree. If there is equality,
+ * choose any vertex of maximal degree in the uncolored subgraph.
+ * -# Color the chosen vertex with the least possible (lowest numbered) color.
+ * -# If all the vertices are colored, stop. Otherwise, return to 3.
+ */
+ template <typename Iterator>
+ std::vector<std::vector<Iterator> > make_dsatur_coloring(std::vector<Iterator> &partition,
+ std::vector<types::global_dof_index> (*get_conflict_indices)(Iterator &))
+ {
+ std::vector<std::vector<Iterator> > partition_coloring;
+ // Number of zones composing the partitioning.
+ const unsigned int partition_size(partition.size());
+ std::vector<unsigned int> sorted_vertices(partition_size);
+ std::vector<unsigned int> degrees(partition_size);
+ std::vector<std::vector<types::global_dof_index> > conflict_indices(partition_size);
+ std::vector<std::vector<unsigned int> > graph(partition_size);
+
+ // Get the conflict indices associated to each iterator. The conflict_indices have to be sorted so
+ // set_intersection can be used later.
+ for (unsigned int i=0; i<partition_size; ++i)
+ {
+ conflict_indices[i] = (*get_conflict_indices)(partition[i]);
+ std::sort(conflict_indices[i].begin(),conflict_indices[i].end());
+ }
+
+ // Compute the degree of each vertex of the graph using the
+ // intersection of the conflict indices.
+ std::vector<types::global_dof_index> conflict_indices_intersection;
+ std::vector<types::global_dof_index>::iterator intersection_it;
+ for (unsigned int i=0; i<partition_size; ++i)
+ for (unsigned int j=i+1; j<partition_size; ++j)
+ {
+ conflict_indices_intersection.resize(std::max(conflict_indices[i].size(),
+ conflict_indices[j].size()));
+ intersection_it = std::set_intersection(conflict_indices[i].begin(),
+ conflict_indices[i].end(),conflict_indices[j].begin(),
+ conflict_indices[j].end(),conflict_indices_intersection.begin());
+ // If the two iterators share indices then we increase the degree of the
+ // vertices and create an ''edge'' in the graph.
+ if (intersection_it!=conflict_indices_intersection.begin())
+ {
+ ++degrees[i];
+ ++degrees[j];
+ graph[i].push_back(j);
+ graph[j].push_back(i);
+ }
+ }
+
+ // Sort the vertices by decreasing degree.
+ std::vector<unsigned int>::iterator degrees_it;
+ for (unsigned int i=0; i<partition_size; ++i)
+ {
+ // Find the largest element.
+ degrees_it = std::max_element(degrees.begin(),degrees.end());
+ sorted_vertices[i] = degrees_it-degrees.begin();
+ // Zero the largest element.
+ *degrees_it = 0;
+ }
+
+ // Color the graph.
+ std::vector<boost::unordered_set<unsigned int> > colors_used;
+ for (unsigned int i=0; i<partition_size; ++i)
+ {
+ const unsigned int current_vertex(sorted_vertices[i]);
+ bool new_color(true);
+ // Try to use an existing color, i.e., try to find a color which is not
+ // associated to one of the vertices linked to current_vertex.
+ // Loop over the color.
+ for (unsigned int j=0; j<partition_coloring.size(); ++j)
+ {
+ // Loop on the vertices linked to current_vertex. If one vertex linked
+ // to current_vertex is already using the color j, this color cannot
+ // be used anymore.
+ bool unused_color(true);
+ for (unsigned int k=0; k<graph[current_vertex].size(); ++k)
+ if (colors_used[j].count(graph[current_vertex][k])==1)
+ {
+ unused_color = false;
+ break;
+ }
+ if (unused_color)
+ {
+ partition_coloring[j].push_back(partition[current_vertex]);
+ colors_used[j].insert(current_vertex);
+ new_color = false;
+ break;
+ }
+ }
+ // Add a new color.
+ if (new_color)
+ {
+ partition_coloring.push_back(std::vector<Iterator> (1,
+ partition[current_vertex]));
+ boost::unordered_set<unsigned int> tmp;
+ tmp.insert(current_vertex);
+ colors_used.push_back(tmp);
+ }
+ }
+
+ return partition_coloring;
+ }
+
+
+
+ /**
+ * Given a partition-coloring graph, gather the colors together. All the
+ * colors on even (resp. odd) partition can be executed simultaneously. This
+ * function tries to create colors of similar number of elements.
+ */
+ template <typename Iterator>
+ std::vector<std::vector<Iterator> >
+ gather_colors(std::vector<std::vector<std::vector<Iterator> > > const &partition_coloring)
+ {
+ std::vector<std::vector<Iterator> > coloring;
+
+ // Count the number of iterators in each color.
+ const unsigned int partition_size(partition_coloring.size());
+ std::vector<std::vector<unsigned int> > colors_counter(partition_size);
+ for (unsigned int i=0; i<partition_size; ++i)
+ {
+ const unsigned int n_colors(partition_coloring[i].size());
+ colors_counter[i].resize(n_colors);
+ for (unsigned int j=0; j<n_colors; ++j)
+ colors_counter[i][j] = partition_coloring[i][j].size();
+ }
+
+ // Find the partition with the largest number of colors for the even partition.
+ unsigned int i_color(0);
+ unsigned int max_even_n_colors(0);
+ const unsigned int colors_size(colors_counter.size());
+ for (unsigned int i=0; i<colors_size; i+=2)
+ {
+ if (max_even_n_colors<colors_counter[i].size())
+ {
+ max_even_n_colors = colors_counter[i].size();
+ i_color = i;
+ }
+ }
+ coloring.resize(max_even_n_colors);
+ for (unsigned int j=0; j<colors_counter[i_color].size(); ++j)
+ coloring[j] = partition_coloring[i_color][j];
+
+ for (unsigned int i=0; i<partition_size; i+=2)
+ {
+ if (i!=i_color)
+ {
+ boost::unordered_set<unsigned int> used_k;
+ for (unsigned int j=0; j<colors_counter[i].size(); ++j)
+ {
+ // Find the color in the current partition with the largest number of
+ // iterators.
+ std::vector<unsigned int>::iterator it;
+ it = std::max_element(colors_counter[i].begin(),colors_counter[i].end());
+ unsigned int min_iterators(-1);
+ unsigned int pos(0);
+ // Find the color of coloring with the least number of colors among
+ // the colors that have not been used yet.
+ for (unsigned int k=0; k<max_even_n_colors; ++k)
+ if (used_k.count(k)==0)
+ if (colors_counter[i_color][k]<min_iterators)
+ {
+ min_iterators = colors_counter[i_color][k];
+ pos = k;
+ }
+ colors_counter[i_color][pos] += *it;
+ // Concatenate the current color with the existing coloring.
+ coloring[pos].insert(coloring[pos].end(),
+ partition_coloring[i][it-colors_counter[i].begin()].begin(),
+ partition_coloring[i][it-colors_counter[i].begin()].end());
+ used_k.insert(pos);
+ // Put the number of iterators to the current color to zero.
+ *it = 0;
+ }
+ }
+ }
+
+ // Do the same thing that we did for the even partitions to the odd
+ // partitions
+ unsigned int max_odd_n_colors(0);
+ for (unsigned int i=1; i<partition_size; i+=2)
+ {
+ if (max_odd_n_colors<colors_counter[i].size())
+ {
+ max_odd_n_colors = colors_counter[i].size();
+ i_color = i;
+ }
+ }
+ coloring.resize(max_even_n_colors+max_odd_n_colors);
+ for (unsigned int j=0; j<colors_counter[i_color].size(); ++j)
+ coloring[max_even_n_colors+j] = partition_coloring[i_color][j];
+
+ for (unsigned int i=1; i<partition_size; i+=2)
+ {
+ if (i!=i_color)
+ {
+ boost::unordered_set<unsigned int> used_k;
+ for (unsigned int j=0; j<colors_counter[i].size(); ++j)
+ {
+ // Find the color in the current partition with the largest number of
+ // iterators.
+ std::vector<unsigned int>::iterator it;
+ it = std::max_element(colors_counter[i].begin(),colors_counter[i].end());
+ unsigned int min_iterators(-1);
+ unsigned int pos(0);
+ // Find the color of coloring with the least number of colors among
+ // the colors that have not been used yet.
+ for (unsigned int k=0; k<max_odd_n_colors; ++k)
+ if (used_k.count(k)==0)
+ if (colors_counter[i_color][k]<min_iterators)
+ {
+ min_iterators = colors_counter[i_color][k];
+ pos = k;
+ }
+ colors_counter[i_color][pos] += *it;
+ // Concatenate the current color with the existing coloring.
+ coloring[max_even_n_colors+pos].insert(coloring[max_even_n_colors+pos].end(),
+ partition_coloring[i][it-colors_counter[i].begin()].begin(),
+ partition_coloring[i][it-colors_counter[i].begin()].end());
+ used_k.insert(pos);
+ // Put the number of iterators to the current color to zero.
+ *it = 0;
+ }
+ }
+ }
+
+ return coloring;
+ }
+
+
+
+ /**
+ * This function creates a coloring given two iterators on the DoFHandler
+ * and a function that return the conflict indices given an iterator. When
+ * using continuous finite elements, the conflict_indices can be the dofs
+ * indices.
+ */
+ template <typename Iterator>
+ std::vector<std::vector<Iterator> >
+ make_graph_coloring(Iterator const &begin,typename identity<Iterator>::type const &end,
+ std::vector<types::global_dof_index> (*get_conflict_indices)(Iterator &))
+ {
+ // Create the partitioning.
+ std::vector<std::vector<Iterator> > partitioning = create_partitioning(begin,end,
+ get_conflict_indices);
+
+ // Color the iterators within each partition.
+ const unsigned int partitioning_size(partitioning.size());
+ std::vector<std::vector<std::vector<Iterator> > > partition_coloring(
+ partitioning_size);
+ for (unsigned int i=0; i<partitioning_size; ++i)
+ {
+ // Compute the coloring of the graph using the DSATUR algorithm
+ partition_coloring[i] = make_dsatur_coloring(partitioning[i],get_conflict_indices);
+ }
+
+ // Gather the colors together.
+ std::vector<std::vector<Iterator> > coloring = gather_colors(partition_coloring);
+
+ return coloring;
+ }
+
+} // End graph_coloring namespace
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+//---------------------------- graph_coloring.h ---------------------------
+// end of #ifndef __deal2__graph_coloring_h
+#endif
+//---------------------------- graph_coloring.h ---------------------------
* p<sub>i</sub>(x)p<sub>j</sub>(y)p<sub>k</sub>(z)</i>.
*/
void compute_index (const unsigned int n,
- unsigned int (&index)[dim]) const;
+ unsigned int (&index)[dim>0?dim:1]) const;
private:
/**
}
}
-
+//forward declaration of the data type for periodic face pairs
+namespace GridTools {template <typename CellIterator> struct PeriodicFacePair;}
namespace parallel
{
get_p4est_tree_to_coarse_cell_permutation() const;
+
/**
* Join faces in the p4est forest due to periodic boundary conditions.
*
+ * The vector can be filled by the function
+ * GridTools::collect_periodic_faces.
+ *
+ * @todo At the moment just default orientation is implemented.
+ *
+ * @note Before this function can be used the triangulation has to be
+ * initialized and must not be refined.
+ * Calling this function more than once is possible, but not recommended:
+ * The function destroys and rebuilds the p4est forest each time it is called.
+ */
+ void
+ add_periodicity
+ (const std::vector<GridTools::PeriodicFacePair<cell_iterator> >&);
+
+ /**
+ * Same as the function above, but takes a different argument.
+ *
* The entries in the std::vector should have the form
* std_cxx1x::tuple<cell1, face_no1, cell2, face_no2>.
*
* The vector can be filled by the function
- * DoFTools::identify_periodic_face_pairs.
- *
+ * GridTools::identify_periodic_face_pairs.
+ *
+ * @note This function can only be used if the faces are in
+ * default orientation.
+ *
* @note Before this function can be used the triangulation has to be
* initialized and must not be refined.
* Calling this function more than once is possible, but not recommended:
cell_iterator, unsigned int> >&);
+
private:
/**
* MPI communicator to be
template <class GridClass> class InterGridMap;
template <int dim, int spacedim> class Mapping;
+namespace GridTools {template <typename CellIterator> struct PeriodicFacePair;}
//TODO: map_support_points_to_dofs should generate a multimap, rather than just a map, since several dofs may be located at the same support point
*/
template<typename FaceIterator>
void
- make_periodicity_constraints (const FaceIterator &face_1,
- const typename identity<FaceIterator>::type &face_2,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask(),
- const bool face_orientation = true,
- const bool face_flip = false,
- const bool face_rotation = false);
+ make_periodicity_constraints
+ (const FaceIterator &face_1,
+ const typename identity<FaceIterator>::type &face_2,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask(),
+ const bool face_orientation = true,
+ const bool face_flip = false,
+ const bool face_rotation = false);
+
/**
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask());
+ make_periodicity_constraints
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const int direction,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
/**
*
* @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
*
- * @author Matthias Maier, 2012
+ * @author Daniel Arndt, 2013, Matthias Maier, 2012
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- dealii::Tensor<1,DH::space_dimension> &offset,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask());
+ make_periodicity_constraints
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const int direction,
+ dealii::Tensor<1,DH::space_dimension> &offset,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
/**
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id b_id,
- const int direction,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask());
+ make_periodicity_constraints
+ (const DH &dof_handler,
+ const types::boundary_id b_id,
+ const int direction,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
/**
*/
template<typename DH>
void
- make_periodicity_constraints (const DH &dof_handler,
- const types::boundary_id b_id,
- const int direction,
- dealii::Tensor<1,DH::space_dimension> &offset,
- dealii::ConstraintMatrix &constraint_matrix,
- const ComponentMask &component_mask = ComponentMask());
+ make_periodicity_constraints
+ (const DH &dof_handler,
+ const types::boundary_id b_id,
+ const int direction,
+ dealii::Tensor<1,DH::space_dimension> &offset,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
+
+ /**
+ * Same as above but the periodicity information is given by
+ * @p periodic_faces. This std::vector can be created by
+ * GridTools::collect_periodic_faces.
+ *
+ * @note For DoFHandler objects that are built on a
+ * parallel::distributed::Triangulation object
+ * parallel::distributed::Triangulation::add_periodicity has to be called
+ * before.
+ */
+ template<typename DH>
+ void
+ make_periodicity_constraints
+ (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+ &periodic_faces,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask = ComponentMask());
//@}
#include <deal.II/base/config.h>
#include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/base/polynomial_space.h>
#include <deal.II/fe/fe_poly_face.h>
DEAL_II_NAMESPACE_OPEN
* and two, the polynomials hence correspond to the usual Lagrange polynomials
* on equidistant points.
*
- * This finite element is the trace space of FE_RaviartThomas on the
- * faces and serves in hybridized methods.
+ * This finite element is the trace space of FE_RaviartThomas on the faces and
+ * serves in hybridized methods, e.g. in combination with the FE_DGQ
+ * element. Its use is demonstrated in the step-51 tutorial program.
*
* @note Since these are only finite elements on faces, only
* FEFaceValues and FESubfaceValues will be able to extract reasonable
/**
* Return whether this element implements its hanging node constraints in
* the new way, which has to be used to make elements "hp compatible".
+ */
+ virtual bool hp_constraints_are_implemented () const;
+
+ /**
+ * Return whether this element dominates the one given as argument when they
+ * meet at a common face, whether it is the other way around, whether
+ * neither dominates, or if either could dominate.
*
- * For the FE_Q class the result is always true (independent of the degree
- * of the element), as it implements the complete set of functions necessary
- * for hp capability.
+ * For a definition of domination, see FiniteElementBase::Domination and in
+ * particular the @ref hp_paper "hp paper".
+ */
+ virtual
+ FiniteElementDomination::Domination
+ compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const;
+
+private:
+ /**
+ * Return vector with dofs per vertex, line, quad, hex.
+ */
+ static std::vector<unsigned int> get_dpo_vector (const unsigned int deg);
+};
+
+
+
+
+/**
+ * A finite element, which is a Legendre on each face (i.e., FE_DGP)
+ * and undefined in the interior of the cells. The basis functions on
+ * the faces are from Polynomials::Legendre.
+ *
+ * This element is used in a hybridized method together with the FE_DGP
+ * element for the interior degrees of freedom.
+ *
+ * @note Since these are only finite elements on faces, only
+ * FEFaceValues and FESubfaceValues will be able to extract reasonable
+ * values from any face polynomial. In order to make the use of
+ * FESystem simpler, FEValues objects will not fail using this finite
+ * element space, but all shape function values extracted will equal
+ * to zero.
+ *
+ * @ingroup fe
+ * @author Martin Kronbichler
+ * @date 2013
+ */
+template <int dim, int spacedim=dim>
+class FE_FaceP : public FE_PolyFace<PolynomialSpace<dim-1>, dim, spacedim>
+{
+public:
+ /**
+ * Constructor for complete basis of polynomials of degree <tt>p</tt>. The
+ * shape functions created using this constructor correspond to Legendre
+ * polynomials in each coordinate direction.
+ */
+ FE_FaceP(unsigned int p);
+
+ /**
+ * Clone method.
+ */
+ virtual FiniteElement<dim,spacedim> *clone() const;
+
+ /**
+ * Return a string that uniquely identifies a finite element. This class
+ * returns <tt>FE_FaceP<dim>(degree)</tt> , with <tt>dim</tt> and
+ * <tt>degree</tt> replaced by appropriate values.
+ */
+ virtual std::string get_name () const;
+
+ /**
+ * Return the matrix interpolating from a face of of one element to the face
+ * of the neighboring element. The size of the matrix is then
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+ * element only provides interpolation matrices for elements of the same
+ * type and FE_Nothing. For all other elements, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+ FullMatrix<double> &matrix) const;
+
+ /**
+ * Return the matrix interpolating from a face of of one element to the face
+ * of the neighboring element. The size of the matrix is then
+ * <tt>source.dofs_per_face</tt> times <tt>this->dofs_per_face</tt>. This
+ * element only provides interpolation matrices for elements of the same
+ * type and FE_Nothing. For all other elements, an exception of type
+ * FiniteElement<dim,spacedim>::ExcInterpolationNotImplemented is thrown.
+ */
+ virtual void
+ get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &source,
+ const unsigned int subface,
+ FullMatrix<double> &matrix) const;
+
+ /**
+ * Check for non-zero values on a face.
+ *
+ * This function returns @p true, if the shape function @p shape_index has
+ * non-zero values on the face @p face_index.
+ *
+ * Implementation of the interface in FiniteElement
+ */
+ virtual bool has_support_on_face (const unsigned int shape_index,
+ const unsigned int face_index) const;
+
+ /**
+ * Return whether this element implements its hanging node constraints in
+ * the new way, which has to be used to make elements "hp compatible".
*/
virtual bool hp_constraints_are_implemented () const;
static std::vector<unsigned int> get_dpo_vector (const unsigned int deg);
};
+
DEAL_II_NAMESPACE_CLOSE
#endif
* a hanging node located on a face or an
* edge of it.
*
+ * @note If the point requested does not lie in any of the cells of
+ * the mesh given, then this function throws an exception of type
+ * GridTools::ExcPointNotFound. You can catch this exception and
+ * decide what to do in that case.
+ *
* @note It isn't entirely clear at this time whether the function
* does the right thing with anisotropically refined meshes. It needs
* to be checked for this case.
/**
- * Find and return an iterator to
- * the active cell that surrounds
- * a given point @p ref. The
- * type of the first parameter
- * may be either
- * Triangulation,
- * DoFHandler, or
- * MGDoFHandler, i.e. we
- * can find the cell around a
- * point for iterators into each
- * of these classes.
- *
- * This is solely a wrapper function
- * for the @p interpolate function
- * given below,
- * providing backward compatibility.
- * A Q1 mapping is used for the
- * boundary, and the iterator to
- * the cell in which the point
- * resides is returned.
- *
- * It is recommended to use the
- * other version of this function,
- * as it simultaneously delivers the
- * local coordinate of the given point
- * without additional computational cost.
+ * Find and return an iterator to the active cell that surrounds a
+ * given point @p ref. The type of the first parameter may be either
+ * Triangulation, or one of the DoF handler classes, i.e. we can find the
+ * cell around a point for iterators into each of these classes.
+ *
+ * This is solely a wrapper function for the function of same name
+ * below. A Q1 mapping is used for the boundary, and the iterator
+ * to the cell in which the point resides is returned.
+ *
+ * It is recommended to use the other version of this function, as
+ * it simultaneously delivers the local coordinate of the given
+ * point without additional computational cost.
+ *
+ * @note If the point requested does not lie in any of the cells of
+ * the mesh given, then this function throws an exception of type
+ * GridTools::ExcPointNotFound. You can catch this exception and
+ * decide what to do in that case.
*
* @note When applied to a triangulation or DoF handler object based
* on a parallel::distributed::Triangulation object, the cell
const Point<spacedim> &p);
/**
- * Find and return an iterator to
- * the active cell that surrounds
- * a given point @p p. The
- * type of the first parameter
- * may be either
- * Triangulation,
- * DoFHandler, hp::DoFHandler, or
- * MGDoFHandler, i.e., we
- * can find the cell around a
- * point for iterators into each
- * of these classes.
- *
- * The algorithm used in this
- * function proceeds by first
- * looking for vertex located
- * closest to the given point, see
- * find_closest_vertex(). Secondly,
- * all adjacent cells to this point
- * are found in the mesh, see
- * find_cells_adjacent_to_vertex().
- * Lastly, for each of these cells,
- * it is tested whether the point is
- * inside. This check is performed
- * using arbitrary boundary mappings.
- * Still, it is possible that due
- * to roundoff errors, the point
- * cannot be located exactly inside
- * the unit cell. In this case,
- * even points at a very small
- * distance outside the unit cell
+ * Find and return an iterator to the active cell that surrounds a
+ * given point @p p. The type of the first parameter may be either
+ * Triangulation, DoFHandler, hp::DoFHandler, or MGDoFHandler, i.e.,
+ * we can find the cell around a point for iterators into each of
+ * these classes.
+ *
+ * The algorithm used in this function proceeds by first looking for
+ * vertex located closest to the given point, see
+ * find_closest_vertex(). Secondly, all adjacent cells to this point
+ * are found in the mesh, see find_cells_adjacent_to_vertex().
+ * Lastly, for each of these cells, it is tested whether the point
+ * is inside. This check is performed using arbitrary boundary
+ * mappings. Still, it is possible that due to roundoff errors, the
+ * point cannot be located exactly inside the unit cell. In this
+ * case, even points at a very small distance outside the unit cell
* are allowed.
*
- * If a point lies on the
- * boundary of two or more cells,
- * then the algorithm tries to identify
- * the cell that is of highest
+ * If a point lies on the boundary of two or more cells, then the
+ * algorithm tries to identify the cell that is of highest
* refinement level.
*
- * The function returns an
- * iterator to the cell, as well
- * as the local position of the
- * point inside the unit
- * cell. This local position
- * might be located slightly
- * outside an actual unit cell,
- * due to numerical roundoff.
- * Therefore, the point returned
- * by this function should
- * be projected onto the unit cell,
- * using GeometryInfo::project_to_unit_cell.
- * This is not automatically performed
- * by the algorithm.
+ * The function returns an iterator to the cell, as well as the
+ * local position of the point inside the unit cell. This local
+ * position might be located slightly outside an actual unit cell,
+ * due to numerical roundoff. Therefore, the point returned by this
+ * function should be projected onto the unit cell, using
+ * GeometryInfo::project_to_unit_cell. This is not automatically
+ * performed by the algorithm.
+ *
+ * @note If the point requested does not lie in any of the cells of
+ * the mesh given, then this function throws an exception of type
+ * GridTools::ExcPointNotFound. You can catch this exception and
+ * decide what to do in that case.
*
* @note When applied to a triangulation or DoF handler object based
* on a parallel::distributed::Triangulation object, the cell
const Point<spacedim> &p);
/**
- * A version of the previous function
- * where we use that mapping on a given
- * cell that corresponds to the active
- * finite element index of that
- * cell. This is obviously only useful
- * for hp problems, since the active
- * finite element index for all other DoF
- * handlers is always zero.
+ * A version of the previous function where we use that mapping on a
+ * given cell that corresponds to the active finite element index of
+ * that cell. This is obviously only useful for hp problems, since
+ * the active finite element index for all other DoF handlers is
+ * always zero.
+ *
+ * @note If the point requested does not lie in any of the cells of
+ * the mesh given, then this function throws an exception of type
+ * GridTools::ExcPointNotFound. You can catch this exception and
+ * decide what to do in that case.
*
* @note When applied to a triangulation or DoF handler object based
* on a parallel::distributed::Triangulation object, the cell
const std::set<types::boundary_id> &boundary_ids
= std::set<types::boundary_id>());
-
+ /**
+ * Data type that provides all the information that is needed
+ * to create periodicity constraints and a periodic p4est forest
+ * with respect to two periodic cell faces
+ */
+ template<typename CellIterator>
+ struct PeriodicFacePair
+ {
+ CellIterator cell[2];
+ unsigned int face_idx[2];
+ std::bitset<3> orientation;
+ };
/**
* An orthogonal equality test for faces.
/**
- * This function loops over all faces from @p begin to @p end and
- * collects a set of periodic face pairs for @p direction:
+ * This function will collect periodic face pairs on the highest (i.e.
+ * coarsest) mesh level.
*
* Define a 'first' boundary as all boundary faces having boundary_id
* @p b_id1 and a 'second' boundary consisting of all faces belonging
* boundary with faces belonging to the second boundary with the help
* of orthogonal_equality().
*
- * The bitset that is returned together with the second face encodes the
+ * The bitset that is returned inside of PeriodicFacePair encodes the
* _relative_ orientation of the first face with respect to the second
* face, see the documentation of orthogonal_equality for further details.
*
+ * The @p direction refers to the space direction in which periodicity
+ * is enforced.
+ *
* The @p offset is a vector tangential to the faces that is added to the
* location of vertices of the 'first' boundary when attempting to match
* them to the corresponding vertices of the 'second' boundary. This can
* be used to implement conditions such as $u(0,y)=u(1,y+1)$.
*
- * @author Matthias Maier, 2012
- */
- template<typename FaceIterator>
- std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
- collect_periodic_face_pairs (const FaceIterator &begin,
- const typename identity<FaceIterator>::type &end,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset);
-
-
- /**
- * Same function as above, but accepts a Triangulation or DoFHandler
- * object @p container (a container is a collection of objects, here a
- * collection of cells) instead of an explicit face iterator range.
+ * @note The created std::vector can be used in
+ * DoFTools::make_periodicity_constraints and in
+ * parallel::distributed::Triangulation::add_periodicity to enforce
+ * periodicity algebraically.
*
- * This function will collect periodic face pairs on the highest (i.e.
- * coarsest) mesh level.
- *
- * @author Matthias Maier, 2012
+ * @author Daniel Arndt, 2013
*/
template<typename DH>
- std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
- collect_periodic_face_pairs (const DH &container,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- const dealii::Tensor<1,DH::space_dimension> &offset);
-
+ std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ collect_periodic_faces
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const unsigned int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset);
/**
* This compatibility version of collect_periodic_face_pairs only works
* meshes with cells not in @ref GlossFaceOrientation
* "standard orientation".
*
+ * @author Daniel Arndt, 2013
+ */
+ template<typename DH>
+ std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ collect_periodic_faces
+ (const DH &dof_handler,
+ const types::boundary_id b_id,
+ const unsigned int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset);
+
+ /**
+ * This function does the same as collect_periodic_faces but returns a
+ * different data type.
+ *
* @author Matthias Maier, 2012
+ *
+ * @note The returned data type is not compatible with
+ * DoFTools::make_periodicity_constraints
+ *
+ * @deprecated
+ */
+ template<typename DH>
+ std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
+ collect_periodic_face_pairs
+ (const DH &container,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset) DEAL_II_DEPRECATED;
+
+ /**
+ * This compatibility version of collect_periodic_face_pairs only works
+ * on grids with cells in @ref GlossFaceOrientation "standard orientation".
+ *
+ * @author Matthias Maier, 2012
+ *
+ * @note The returned data type is not compatible with
+ * DoFTools::make_periodicity_constraints
+ *
+ * @deprecated
*/
template<typename DH>
std::map<typename DH::face_iterator, typename DH::face_iterator>
- collect_periodic_face_pairs (const DH &dof_handler, /*TODO: Name*/
- const types::boundary_id b_id,
- const int direction,
- const dealii::Tensor<1,DH::space_dimension> &offset);
+ collect_periodic_face_pairs
+ (const DH &dof_handler, /*TODO: Name*/
+ const types::boundary_id b_id,
+ int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset) DEAL_II_DEPRECATED;
/**
- * Add periodicity information to the @p periodicity_vector for the
- * boundaries with boundary id @p b_id1 and @p b_id2 in cartesian
- * direction @p direction.
+ * This version loops over all faces from @p begin to @p end
+ * instead of accepting a DoFHandler or a Triangulation.
*
- * This function tries to match all faces belonging to the first
- * boundary with faces belonging to the second boundary
- * by comparing the center of the cell faces. To find the correct
- * corresponding faces, the direction argument indicates in which
- * cartesian direction periodicity should be set.
- * ((0,1,2) -> (x,y,z)-direction)
+ * @author Matthias Maier, 2012
+ *
+ * @note This function cannot produce the return as the other
+ * collect_periodic_faces functions.
+ *
+ * @deprecated
+ */
+ template<typename FaceIterator>
+ std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
+ collect_periodic_face_pairs
+ (const FaceIterator &begin,
+ const typename identity<FaceIterator>::type &end,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ DEAL_II_DEPRECATED;
+
+ /**
+ * This function does the same as collect_periodic_faces but returns a
+ * different data type.
*
- * The output of this function can be used in
+ * @author Daniel Arndt, 2013
+ *
+ * @note The returned data type is not compatible with
+ * DoFTools::make_periodicity_constraints, but with a version of
* parallel::distributed::Triangulation::add_periodicity
+ *
+ * @note Use collect_periodic_faces instead.
+ *
+ * @deprecated
*/
template<typename DH>
void
identify_periodic_face_pairs
- (const DH &dof_handler,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const unsigned int direction,
- std::vector<std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
- typename DH::cell_iterator, unsigned int> >
- &periodicity_vector);
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const unsigned int direction,
+ std::vector<std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+ typename DH::cell_iterator, unsigned int> >
+ &periodicity_vector) DEAL_II_DEPRECATED;
/**
// make sure we are on an
// object for which DoFs have
// been allocated at all
- if (dof_offsets[obj_index] == numbers::invalid_dof_index)
+ if (dof_offsets[obj_index] == numbers::invalid_unsigned_int)
return 0;
Assert (structdim<dim, ExcMessage ("This object can not be used for cells."));
// make sure we are on an
// object for which DoFs have
// been allocated at all
- Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ Assert (dof_offsets[obj_index] != (offset_type)(-1),
ExcMessage ("You are trying to access degree of freedom "
"information for an object on which no such "
"information is available"));
// make sure we are on an
// object for which DoFs have
// been allocated at all
- Assert (dof_offsets[obj_index] != numbers::invalid_dof_index,
+ Assert (dof_offsets[obj_index] != (offset_type)(-1),
ExcMessage ("You are trying to access degree of freedom "
"information for an object on which no such "
"information is available"));
const MDt &Dt,
const MC &C,
VectorMemory<BlockVector<double> > &mem,
- const std::vector<types::global_dof_index> &signature)
+ const std::vector<unsigned int> &signature)
: Ainv(&Ainv), B(&B), Dt(&Dt), C(&C),
mem(mem),
signature(signature),
*/
std::vector<double> vertices;
std::vector<unsigned int> cells;
- unsigned int cell_offset;
+ unsigned int node_offset;
};
vertices.resize(local_points_cell_count[0]*dim);
cells.resize(local_points_cell_count[1]*entries_per_cell);
- cell_offset = global_points_cell_offsets[1]*entries_per_cell;
+ node_offset = global_points_cell_offsets[0];
}
template<int dim>
unsigned int d3)
{
unsigned int base_entry = index * GeometryInfo<dim>::vertices_per_cell;
- cells[base_entry+0] = cell_offset+start;
- cells[base_entry+1] = cell_offset+start+d1;
+ cells[base_entry+0] = node_offset+start;
+ cells[base_entry+1] = node_offset+start+d1;
if (dim>=2)
{
- cells[base_entry+2] = cell_offset+start+d2+d1;
- cells[base_entry+3] = cell_offset+start+d2;
+ cells[base_entry+2] = node_offset+start+d2+d1;
+ cells[base_entry+3] = node_offset+start+d2;
if (dim>=3)
{
- cells[base_entry+4] = cell_offset+start+d3;
- cells[base_entry+5] = cell_offset+start+d3+d1;
- cells[base_entry+6] = cell_offset+start+d3+d2+d1;
- cells[base_entry+7] = cell_offset+start+d3+d2;
+ cells[base_entry+4] = node_offset+start+d3;
+ cells[base_entry+5] = node_offset+start+d3+d1;
+ cells[base_entry+6] = node_offset+start+d3+d2+d1;
+ cells[base_entry+7] = node_offset+start+d3+d2;
}
}
}
while (n_th_vector < vector_data_ranges.size() && std_cxx1x::get<0>(vector_data_ranges[n_th_vector]) < data_set) n_th_vector++;
// Determine whether the data is multiple dimensions or one
- if (std_cxx1x::get<0>(vector_data_ranges[n_th_vector]) == data_set)
+ if (n_th_vector < vector_data_ranges.size() && std_cxx1x::get<0>(vector_data_ranges[n_th_vector]) == data_set)
{
// Multiple dimensions
pt_data_vector_dim = std_cxx1x::get<1>(vector_data_ranges[n_th_vector]) - std_cxx1x::get<0>(vector_data_ranges[n_th_vector])+1;
while (n_th_vector < vector_data_ranges.size() && std_cxx1x::get<0>(vector_data_ranges[n_th_vector]) < data_set) n_th_vector++;
// Determine the dimension of this data
- if (std_cxx1x::get<0>(vector_data_ranges[n_th_vector]) == data_set)
+ if (n_th_vector < vector_data_ranges.size() && std_cxx1x::get<0>(vector_data_ranges[n_th_vector]) == data_set)
{
// Multiple dimensions
pt_data_vector_dim = std_cxx1x::get<1>(vector_data_ranges[n_th_vector]) - std_cxx1x::get<0>(vector_data_ranges[n_th_vector])+1;
ExcMessage ("You can only create a single object of this class "
"in a program since it initializes the MPI system."));
-#ifdef DEAL_II_WITH_PETSC
-# ifdef DEAL_II_WITH_SLEPC
- // Initialise SLEPc (with PETSc):
- SlepcInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
-# else
- // or just initialise PETSc alone:
- PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
-# endif
-#else
-# ifdef DEAL_II_WITH_MPI
+
+
+#ifdef DEAL_II_WITH_MPI
// if we have PETSc, we will initialize it and let it handle MPI.
// Otherwise, we will do it.
int MPI_has_been_started = 0;
AssertThrow (MPI_has_been_started == 0,
ExcMessage ("MPI error. You can only start MPI once!"));
- int mpi_err;
- mpi_err = MPI_Init (&argc, &argv);
+ int mpi_err, provided;
+ // this works likempi_err = MPI_Init (&argc, &argv); but tells MPI that
+ // we might use several threads but never call two MPI functions at the
+ // same time. For an explanation see on why we do this see
+ // http://www.open-mpi.org/community/lists/users/2010/03/12244.php
+ int wanted = MPI_THREAD_SERIALIZED;
+ mpi_err = MPI_Init_thread(&argc, &argv, wanted, &provided);
AssertThrow (mpi_err == 0,
ExcMessage ("MPI could not be initialized."));
-# else
+
+ // disable for now because at least some implementations always return MPI_THREAD_SINGLE.
+ //Assert(max_num_threads==1 || provided != MPI_THREAD_SINGLE,
+ // ExcMessage("MPI reports that we are not allowed to use multiple threads."));
+#else
// make sure the compiler doesn't warn
// about these variables
(void)argc;
(void)argv;
(void)owns_mpi;
+#endif
+
+ // we are allowed to call MPI_Init ourselves and PETScInitialize will
+ // detect this. This allows us to use MPI_Init_thread instead.
+#ifdef DEAL_II_WITH_PETSC
+# ifdef DEAL_II_WITH_SLEPC
+ // Initialize SLEPc (with PETSc):
+ SlepcInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
+# else
+ // or just initialize PETSc alone:
+ PetscInitialize(&argc, &argv, PETSC_NULL, PETSC_NULL);
# endif
#endif
#endif
- // only MPI_Finalize if we are running with MPI and we are not using PETSc
- // (otherwise we called it above already):
+ // only MPI_Finalize if we are running with MPI. We also need to do this
+ // when running PETSc, because we initialize MPI ourselves before calling
+ // PetscInitialize
#ifdef DEAL_II_WITH_MPI
-#ifndef DEAL_II_WITH_PETSC
- int mpi_err = 0;
+ int mpi_err = 0;
int MPI_has_been_started = 0;
MPI_Initialized(&MPI_has_been_started);
AssertThrow (mpi_err == 0,
ExcMessage ("An error occurred while calling MPI_Finalize()"));
-#endif
#endif
}
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/distributed/tria.h>
#include <deal.II/grid/grid_tools.h>
+#include <deal.II/distributed/tria.h>
#ifdef DEAL_II_WITH_P4EST
# include <p4est_bits.h>
return mpi_communicator;
}
-
- template <int dim, int spacedim>
+ template<int dim, int spacedim>
void
Triangulation<dim,spacedim>::add_periodicity
- (const std::vector<std_cxx1x::tuple<cell_iterator, unsigned int,
- cell_iterator, unsigned int> >&
- periodicity_vector)
+ (const std::vector<GridTools::PeriodicFacePair<cell_iterator> >&
+ periodicity_vector)
{
-#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+ #if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
Assert (triangulation_has_content == true,
ExcMessage ("The triangulation is empty!"));
Assert (this->n_levels() == 1,
ExcMessage ("The triangulation is refined!"));
-
- typename std::vector
- <typename std_cxx1x::tuple
- <cell_iterator, unsigned int,
- cell_iterator, unsigned int> >::const_iterator periodic_tuple;
- periodic_tuple = periodicity_vector.begin();
-
- typename std::vector
- <typename std_cxx1x::tuple
- <cell_iterator, unsigned int,
- cell_iterator, unsigned int> >::const_iterator periodic_end;
+
+ typedef std::vector<GridTools::PeriodicFacePair<cell_iterator> >
+ FaceVector;
+ typename FaceVector::const_iterator it, periodic_end;
+ it = periodicity_vector.begin();
periodic_end = periodicity_vector.end();
-
- for (; periodic_tuple<periodic_end; ++periodic_tuple)
+
+ for (; it<periodic_end; ++it)
{
- const cell_iterator first_cell=std_cxx1x::get<0>(*periodic_tuple);
- const cell_iterator second_cell=std_cxx1x::get<2>(*periodic_tuple);
- const unsigned int face_right=std_cxx1x::get<3>(*periodic_tuple);
- const unsigned int face_left=std_cxx1x::get<1>(*periodic_tuple);
-
+ const cell_iterator first_cell = it->cell[0];
+ const cell_iterator second_cell = it->cell[1];
+ const unsigned int face_left = it->face_idx[0];
+ const unsigned int face_right = it->face_idx[1];
+
//respective cells of the matching faces in p4est
const unsigned int tree_left
= coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
= coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
second_cell)];
+ //TODO Add support for non default orientation.
+ Assert(it->orientation == 1,
+ ExcMessage("Found a face match with non standard orientation. "
+ "This function is only suitable for meshes with "
+ "cells in default orientation"));
+
dealii::internal::p4est::functions<dim>::
connectivity_join_faces (connectivity,
tree_left,
face_left,
face_right,
/* orientation */ 0);
+
/* The orientation parameter above describes the difference between
* the cell on the left and the cell on the right would number of the
* corners of the face. In the periodic domains most users will want,
* date.
*/
}
-
-
+
+
Assert(dealii::internal::p4est::functions<dim>::connectivity_is_valid
- (connectivity) == 1,
- ExcInternalError());
-
- // now create a forest out of the connectivity data structure
+ (connectivity) == 1, ExcInternalError());
+
+ // now create a forest out of the connectivity data structure
dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
parallel_forest
= dealii::internal::p4est::functions<dim>::
new_forest (mpi_communicator,
connectivity,
- /* minimum initial number of quadrants per tree */ 0,
- /* minimum level of upfront refinement */ 0,
- /* use uniform upfront refinement */ 1,
- /* user_data_size = */ 0,
- /* user_data_constructor = */ NULL,
- /* user_pointer */ this);
+ /* minimum initial number of quadrants per tree */ 0,
+ /* minimum level of upfront refinement */ 0,
+ /* use uniform upfront refinement */ 1,
+ /* user_data_size = */ 0,
+ /* user_data_constructor = */ NULL,
+ /* user_pointer */ this);
+
+ #else
+ Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
+ #endif
+ }
+
+
+ template <int dim, int spacedim>
+ void
+ Triangulation<dim,spacedim>::add_periodicity
+ (const std::vector<std_cxx1x::tuple<cell_iterator, unsigned int,
+ cell_iterator, unsigned int> >&
+ periodicity_vector)
+ {
+#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+ typedef std::vector<std_cxx1x::tuple <cell_iterator, unsigned int,
+ cell_iterator, unsigned int> >
+ FaceVector;
+
+ typename FaceVector::const_iterator it, end_periodic;
+ it = periodicity_vector.begin();
+ end_periodic = periodicity_vector.end();
+
+ std::vector<GridTools::PeriodicFacePair<cell_iterator> > periodic_faces;
+
+ for(; it!=end_periodic; ++it)
+ {
+ const cell_iterator cell1 = std_cxx1x::get<0> (*it);
+ const cell_iterator cell2 = std_cxx1x::get<2> (*it);
+ const unsigned int face_idx1 = std_cxx1x::get<1> (*it);
+ const unsigned int face_idx2 = std_cxx1x::get<3> (*it);
+ const GridTools::PeriodicFacePair<cell_iterator> matched_face
+ = {{cell1, cell2},{face_idx1, face_idx2}, 1};
+ periodic_faces.push_back(matched_face);
+ }
+ add_periodicity(periodic_faces);
#else
Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
#endif
}
-
+
template<typename DH>
void
make_periodicity_constraints (const DH &dof_handler,
ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
"different to denote different boundaries."));
- typedef typename DH::face_iterator FaceIterator;
- typedef std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > > FaceMap;
+ typedef std::vector<GridTools::PeriodicFacePair
+ <typename DH::cell_iterator> > FaceVector;
// Collect matching periodic cells on the coarsest level:
- FaceMap matched_cells =
- GridTools::collect_periodic_face_pairs(dof_handler,
- b_id1, b_id2,
- direction, offset);
-
- // And apply the low level make_periodicity_constraints function to
- // every matching pair:
- for (typename FaceMap::iterator it = matched_cells.begin();
- it != matched_cells.end(); ++it)
- {
- typedef typename DH::face_iterator FaceIterator;
- const FaceIterator &face_1 = it->first;
- const FaceIterator &face_2 = it->second.first;
- const std::bitset<3> &orientation = it->second.second;
-
- Assert(face_1->at_boundary() && face_2->at_boundary(),
- ExcInternalError());
-
- Assert (face_1->boundary_indicator() == b_id1 &&
- face_2->boundary_indicator() == b_id2,
- ExcInternalError());
-
- Assert (face_1 != face_2,
- ExcInternalError());
+ FaceVector matched_faces =
+ GridTools::collect_periodic_faces(dof_handler, b_id1, b_id2,
+ direction, offset);
- make_periodicity_constraints(face_1,
- face_2,
- constraint_matrix,
- component_mask,
- orientation[0],
- orientation[1],
- orientation[2]);
- }
+ make_periodicity_constraints<DH>
+ (matched_faces, constraint_matrix, component_mask);
}
Assert(dim == space_dim,
ExcNotImplemented());
- typedef typename DH::face_iterator FaceIterator;
- typedef std::map<FaceIterator, FaceIterator> FaceMap;
+ typedef std::vector<GridTools::PeriodicFacePair
+ <typename DH::cell_iterator> > FaceVector;
// Collect matching periodic cells on the coarsest level:
- FaceMap matched_cells =
- GridTools::collect_periodic_face_pairs(dof_handler,
- b_id,
- direction, offset);
+ const FaceVector matched_faces =
+ GridTools::collect_periodic_faces(dof_handler, b_id,
+ direction, offset);
+
+ make_periodicity_constraints<DH>
+ (matched_faces, constraint_matrix, component_mask);
+ }
+
+
+
+ template<typename DH>
+ void
+ make_periodicity_constraints
+ (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+ &periodic_faces,
+ dealii::ConstraintMatrix &constraint_matrix,
+ const ComponentMask &component_mask)
+ {
+ typedef std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+ FaceVector;
+ typename FaceVector::const_iterator it, end_periodic;
+ it = periodic_faces.begin();
+ end_periodic = periodic_faces.end();
+
// And apply the low level make_periodicity_constraints function to
// every matching pair:
- for (typename FaceMap::iterator it = matched_cells.begin();
- it != matched_cells.end(); ++it)
- {
- typedef typename DH::face_iterator FaceIterator;
- const FaceIterator &face_1 = it->first;
- const FaceIterator &face_2 = it->second;
-
- Assert(face_1->at_boundary() && face_2->at_boundary(),
- ExcInternalError());
+ for(; it!=end_periodic; ++it)
+ {
+ typedef typename DH::face_iterator FaceIterator;
+ const FaceIterator face_1 = it->cell[0]->face(it->face_idx[0]);
+ const FaceIterator face_2 = it->cell[1]->face(it->face_idx[1]);
- Assert (face_1->boundary_indicator() == b_id &&
- face_2->boundary_indicator() == b_id,
- ExcInternalError());
+ Assert(face_1->at_boundary() && face_2->at_boundary(),
+ ExcInternalError());
- Assert (face_1 != face_2,
- ExcInternalError());
+ Assert (face_1 != face_2,
+ ExcInternalError());
- make_periodicity_constraints(face_1,
- face_2,
- constraint_matrix,
- component_mask
- /* standard orientation */);
- }
+ make_periodicity_constraints(face_1,
+ face_2,
+ constraint_matrix,
+ component_mask,
+ it->orientation[0],
+ it->orientation[1],
+ it->orientation[2]);
+ }
}
dealii::Tensor<1,DH::space_dimension> &,
dealii::ConstraintMatrix &,
const ComponentMask &);
+
+ template
+ void
+ DoFTools::make_periodicity_constraints<DH>
+ (const std::vector<GridTools::PeriodicFacePair<DH::cell_iterator> > &,
+ dealii::ConstraintMatrix &,
+ const ComponentMask &);
+
#endif
}
std::vector<ComponentMask>(FiniteElementData<dim>(
get_dpo_vector(degree), 1, degree).dofs_per_cell, std::vector<bool>(1,true)))
{
- // Reinit the vectors of
- // restriction and prolongation
- // matrices to the right sizes
+ // Reinit the vectors of restriction and prolongation matrices to the right
+ // sizes
this->reinit_restriction_and_prolongation_matrices();
// Fill prolongation matrices with embedding operators
if (dim == spacedim)
std::string
FE_DGP<dim,spacedim>::get_name () const
{
- // note that the
- // FETools::get_fe_from_name
- // function depends on the
- // particular format of the string
- // this function returns, so they
- // have to be kept in synch
+ // note that the FETools::get_fe_from_name function depends on the
+ // particular format of the string this function returns, so they have to be
+ // kept in synch
std::ostringstream namebuf;
namebuf << "FE_DGP<" << dim << ">(" << this->degree << ")";
get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
FullMatrix<double> &interpolation_matrix) const
{
- // this is only implemented, if the source
- // FE is also a DGP element. in that case,
- // both elements have no dofs on their
- // faces and the face interpolation matrix
- // is necessarily empty -- i.e. there isn't
- // much we need to do here.
+ // this is only implemented, if the source FE is also a DGP element. in that
+ // case, both elements have no dofs on their faces and the face
+ // interpolation matrix is necessarily empty -- i.e. there isn't much we
+ // need to do here.
typedef FiniteElement<dim,spacedim> FE;
typedef FE_DGP<dim,spacedim> FEDGP;
AssertThrow ((x_source_fe.get_name().find ("FE_DGP<") == 0)
const unsigned int ,
FullMatrix<double> &interpolation_matrix) const
{
- // this is only implemented, if the source
- // FE is also a DGP element. in that case,
- // both elements have no dofs on their
- // faces and the face interpolation matrix
- // is necessarily empty -- i.e. there isn't
- // much we need to do here.
+ // this is only implemented, if the source FE is also a DGP element. in that
+ // case, both elements have no dofs on their faces and the face
+ // interpolation matrix is necessarily empty -- i.e. there isn't much we
+ // need to do here.
typedef FiniteElement<dim,spacedim> FE;
typedef FE_DGP<dim,spacedim> FEDGP;
AssertThrow ((x_source_fe.get_name().find ("FE_DGP<") == 0)
FE_DGP<dim,spacedim>::
hp_vertex_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
{
- // there are no such constraints for DGP
- // elements at all
+ // there are no such constraints for DGP elements at all
if (dynamic_cast<const FE_DGP<dim,spacedim>*>(&fe_other) != 0)
return
std::vector<std::pair<unsigned int, unsigned int> > ();
FE_DGP<dim,spacedim>::
hp_line_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
{
- // there are no such constraints for DGP
- // elements at all
+ // there are no such constraints for DGP elements at all
if (dynamic_cast<const FE_DGP<dim,spacedim>*>(&fe_other) != 0)
return
std::vector<std::pair<unsigned int, unsigned int> > ();
FE_DGP<dim,spacedim>::
hp_quad_dof_identities (const FiniteElement<dim,spacedim> &fe_other) const
{
- // there are no such constraints for DGP
- // elements at all
+ // there are no such constraints for DGP elements at all
if (dynamic_cast<const FE_DGP<dim,spacedim>*>(&fe_other) != 0)
return
std::vector<std::pair<unsigned int, unsigned int> > ();
FiniteElementDomination::Domination
FE_DGP<dim,spacedim>::compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
{
- // check whether both are discontinuous
- // elements, see the description of
+ // check whether both are discontinuous elements, see the description of
// FiniteElementDomination::Domination
if (dynamic_cast<const FE_DGP<dim,spacedim>*>(&fe_other) != 0)
return FiniteElementDomination::no_requirements;
FE_DGP<dim,spacedim>::has_support_on_face (const unsigned int,
const unsigned int) const
{
- // all shape functions have support on all
- // faces
+ // all shape functions have support on all faces
return true;
}
#include <deal.II/fe/fe_face.h>
#include <deal.II/fe/fe_poly_face.templates.h>
#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_tools.h>
#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/lac/householder.h>
#include <sstream>
DEAL_II_NAMESPACE_OPEN
AssertDimension (k, this->unit_face_support_points.size());
}
- // initialize unit support points
+ // initialize unit support points (this makes it possible to assign initial
+ // values to FE_FaceQ)
this->unit_support_points.resize(GeometryInfo<dim>::faces_per_cell*
this->unit_face_support_points.size());
const unsigned int n_face_dofs = this->unit_face_support_points.size();
}
+
template <int dim, int spacedim>
FiniteElement<dim,spacedim> *
FE_FaceQ<dim,spacedim>::clone() const
}
+
template <int dim, int spacedim>
std::string
FE_FaceQ<dim,spacedim>::get_name () const
{
- // note that the
- // FETools::get_fe_from_name
- // function depends on the
- // particular format of the string
- // this function returns, so they
- // have to be kept in synch
-
+ // note that the FETools::get_fe_from_name function depends on the
+ // particular format of the string this function returns, so they have to be
+ // kept in synch
std::ostringstream namebuf;
namebuf << "FE_FaceQ<" << dim << ">(" << this->degree << ")";
template <int dim, int spacedim>
void
FE_FaceQ<dim,spacedim>::
-get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source_fe,
FullMatrix<double> &interpolation_matrix) const
{
- // this function is similar to the respective method in FE_Q
-
- // this is only implemented, if the source FE is also a FE_FaceQ element
- AssertThrow ((dynamic_cast<const FE_FaceQ<dim,spacedim> *>(&x_source_fe) != 0),
- (typename FiniteElement<dim,spacedim>::
- ExcInterpolationNotImplemented()));
-
- Assert (interpolation_matrix.n() == this->dofs_per_face,
- ExcDimensionMismatch (interpolation_matrix.n(),
- this->dofs_per_face));
- Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
- ExcDimensionMismatch (interpolation_matrix.m(),
- x_source_fe.dofs_per_face));
-
- // ok, source is a FaceQ element, so we will be able to do the work
- const FE_FaceQ<dim,spacedim> &source_fe
- = dynamic_cast<const FE_FaceQ<dim,spacedim>&>(x_source_fe);
-
- // Make sure that the element for which the DoFs should be constrained is
- // the one with the higher polynomial degree. Actually the procedure will
- // work also if this assertion is not satisfied. But the matrices produced
- // in that case might lead to problems in the hp procedures, which use this
- // method.
- Assert (this->dofs_per_face <= source_fe.dofs_per_face,
- (typename FiniteElement<dim,spacedim>::
- ExcInterpolationNotImplemented ()));
-
- // generate a quadrature with the unit face support points.
- const Quadrature<dim-1> face_quadrature (source_fe.get_unit_face_support_points ());
-
- // Rule of thumb for FP accuracy, that can be expected for a given
- // polynomial degree. This value is used to cut off values close to zero.
- const double eps = 2e-13*(this->degree+1)*(dim-1);
-
- // compute the interpolation matrix by simply taking the value at the
- // support points.
- for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
- {
- const Point<dim-1> &p = face_quadrature.point (i);
-
- for (unsigned int j=0; j<this->dofs_per_face; ++j)
- {
- double matrix_entry = this->poly_space.compute_value (j, p);
-
- // Correct the interpolated value. I.e. if it is close to 1 or 0,
- // make it exactly 1 or 0. Unfortunately, this is required to avoid
- // problems with higher order elements.
- if (std::fabs (matrix_entry - 1.0) < eps)
- matrix_entry = 1.0;
- if (std::fabs (matrix_entry) < eps)
- matrix_entry = 0.0;
-
- interpolation_matrix(i,j) = matrix_entry;
- }
- }
-
- // make sure that the row sum of each of the matrices is 1 at this
- // point. this must be so since the shape functions sum up to 1
- for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
- {
- double sum = 0.;
-
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- sum += interpolation_matrix(j,i);
-
- Assert (std::fabs(sum-1) < eps, ExcInternalError());
- }
+ get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int,
+ interpolation_matrix);
}
for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
{
const Point<dim-1> p =
+ subface == numbers::invalid_unsigned_int
+ ?
+ face_quadrature.point(i)
+ :
GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(i),
subface);
return FiniteElementDomination::neither_element_dominates;
}
+
+
+// --------------------------------------- FE_FaceP --------------------------
+
+template <int dim, int spacedim>
+FE_FaceP<dim,spacedim>::FE_FaceP (const unsigned int degree)
+ :
+ FE_PolyFace<PolynomialSpace<dim-1>, dim, spacedim>
+ (PolynomialSpace<dim-1>(Polynomials::Legendre::generate_complete_basis(degree)),
+ FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
+ std::vector<bool>(1,true))
+{}
+
+
+template <int dim, int spacedim>
+FiniteElement<dim,spacedim> *
+FE_FaceP<dim,spacedim>::clone() const
+{
+ return new FE_FaceP<dim,spacedim>(this->degree);
+}
+
+
+template <int dim, int spacedim>
+std::string
+FE_FaceP<dim,spacedim>::get_name () const
+{
+ // note that the FETools::get_fe_from_name function depends on the
+ // particular format of the string this function returns, so they have to be
+ // kept in synch
+ std::ostringstream namebuf;
+ namebuf << "FE_FaceP<" << dim << ">(" << this->degree << ")";
+
+ return namebuf.str();
+}
+
+
+
+template <int dim, int spacedim>
+bool
+FE_FaceP<dim,spacedim>::has_support_on_face (
+ const unsigned int shape_index,
+ const unsigned int face_index) const
+{
+ return (face_index == (shape_index/this->dofs_per_face));
+}
+
+
+
+template <int dim, int spacedim>
+std::vector<unsigned int>
+FE_FaceP<dim,spacedim>::get_dpo_vector (const unsigned int deg)
+{
+ std::vector<unsigned int> dpo(dim+1, 0U);
+ dpo[dim-1] = deg+1;
+ for (unsigned int i=1; i<dim-1; ++i)
+ {
+ dpo[dim-1] *= deg+1+i;
+ dpo[dim-1] /= i+1;
+ }
+ return dpo;
+}
+
+
+
+
+template <int dim, int spacedim>
+bool
+FE_FaceP<dim,spacedim>::hp_constraints_are_implemented () const
+{
+ return true;
+}
+
+
+
+template <int dim, int spacedim>
+FiniteElementDomination::Domination
+FE_FaceP<dim,spacedim>::
+compare_for_face_domination (const FiniteElement<dim,spacedim> &fe_other) const
+{
+ if (const FE_FaceP<dim,spacedim> *fe_q_other
+ = dynamic_cast<const FE_FaceP<dim,spacedim>*>(&fe_other))
+ {
+ if (this->degree < fe_q_other->degree)
+ return FiniteElementDomination::this_element_dominates;
+ else if (this->degree == fe_q_other->degree)
+ return FiniteElementDomination::either_element_can_dominate;
+ else
+ return FiniteElementDomination::other_element_dominates;
+ }
+ else if (dynamic_cast<const FE_Nothing<dim>*>(&fe_other) != 0)
+ {
+ // the FE_Nothing has no degrees of freedom and it is typically used in
+ // a context where we don't require any continuity along the interface
+ return FiniteElementDomination::no_requirements;
+ }
+
+ Assert (false, ExcNotImplemented());
+ return FiniteElementDomination::neither_element_dominates;
+}
+
+
+
+
+template <int dim, int spacedim>
+void
+FE_FaceP<dim,spacedim>::
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source_fe,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int,
+ interpolation_matrix);
+}
+
+
+
+template <int dim, int spacedim>
+void
+FE_FaceP<dim,spacedim>::
+get_subface_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+ const unsigned int subface,
+ FullMatrix<double> &interpolation_matrix) const
+{
+ // this function is similar to the respective method in FE_Q
+
+ Assert (interpolation_matrix.n() == this->dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.n(),
+ this->dofs_per_face));
+ Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
+ ExcDimensionMismatch (interpolation_matrix.m(),
+ x_source_fe.dofs_per_face));
+
+ // see if source is a FaceP element
+ if (const FE_FaceP<dim,spacedim> *source_fe
+ = dynamic_cast<const FE_FaceP<dim,spacedim> *>(&x_source_fe))
+ {
+ // Make sure that the element for which the DoFs should be constrained
+ // is the one with the higher polynomial degree. Actually the procedure
+ // will work also if this assertion is not satisfied. But the matrices
+ // produced in that case might lead to problems in the hp procedures,
+ // which use this method.
+ Assert (this->dofs_per_face <= source_fe->dofs_per_face,
+ (typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented ()));
+
+ // do this as in FETools by solving a least squares problem where we
+ // force the source FE polynomial to be equal the given FE on all
+ // quadrature points
+ const QGauss<dim-1> face_quadrature (source_fe->degree+1);
+
+ // Rule of thumb for FP accuracy, that can be expected for a given
+ // polynomial degree. This value is used to cut off values close to
+ // zero.
+ const double eps = 2e-13*(this->degree+1)*(dim-1);
+
+ FullMatrix<double> mass (face_quadrature.size(), source_fe->dofs_per_face);
+
+ for (unsigned int k = 0; k < face_quadrature.size(); ++k)
+ {
+ const Point<dim-1> p =
+ subface == numbers::invalid_unsigned_int ?
+ face_quadrature.point(k) :
+ GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(k),
+ subface);
+
+ for (unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
+ mass (k , j) = source_fe->poly_space.compute_value(j, p);
+ }
+
+ Householder<double> H(mass);
+ Vector<double> v_in(face_quadrature.size());
+ Vector<double> v_out(source_fe->dofs_per_face);
+
+
+ // compute the interpolation matrix by evaluating on the fine side and
+ // then solving the least squares problem
+ for (unsigned int i=0; i<this->dofs_per_face; ++i)
+ {
+ for (unsigned int k = 0; k < face_quadrature.size(); ++k)
+ {
+ const Point<dim-1> p = numbers::invalid_unsigned_int ?
+ face_quadrature.point(k) :
+ GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(k),
+ subface);
+ v_in(k) = this->poly_space.compute_value(i, p);
+ }
+ const double result = H.least_squares(v_out, v_in);
+ Assert(result < 1e-12, FETools::ExcLeastSquaresError (result));
+
+ for (unsigned int j = 0; j < source_fe->dofs_per_face; ++j)
+ {
+ double matrix_entry = v_out(j);
+
+ // Correct the interpolated value. I.e. if it is close to 1 or 0,
+ // make it exactly 1 or 0. Unfortunately, this is required to avoid
+ // problems with higher order elements.
+ if (std::fabs (matrix_entry - 1.0) < eps)
+ matrix_entry = 1.0;
+ if (std::fabs (matrix_entry) < eps)
+ matrix_entry = 0.0;
+
+ interpolation_matrix(j,i) = matrix_entry;
+ }
+ }
+ }
+ else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != 0)
+ {
+ // nothing to do here, the FE_Nothing has no degrees of freedom anyway
+ }
+ else
+ AssertThrow (false,(typename FiniteElement<dim,spacedim>::
+ ExcInterpolationNotImplemented()));
+}
+
+
// explicit instantiations
#include "fe_face.inst"
{
#if deal_II_dimension > 1
template class FE_PolyFace<TensorProductPolynomials<deal_II_dimension-1> >;
- //template class FE_PolyFace<PolynomialSpace<deal_II_dimension>, deal_II_dimension>;
+ template class FE_PolyFace<PolynomialSpace<deal_II_dimension-1>, deal_II_dimension>;
//template class FE_PolyFace<PolynomialsP<deal_II_dimension>, deal_II_dimension>;
template class FE_FaceQ<deal_II_dimension>;
+ template class FE_FaceP<deal_II_dimension>;
#endif
}
template <class POLY, int dim, int spacedim>
void
FE_Q_Base<POLY,dim,spacedim>::
-get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &x_source_fe,
+get_face_interpolation_matrix (const FiniteElement<dim,spacedim> &source_fe,
FullMatrix<double> &interpolation_matrix) const
{
Assert (dim > 1, ExcImpossibleInDim(1));
-
- // this is only implemented, if the source FE is also a Q element
- AssertThrow ((dynamic_cast<const FE_Q_Base<POLY,dim,spacedim> *>(&x_source_fe) != 0),
- (typename FiniteElement<dim,spacedim>::
- ExcInterpolationNotImplemented()));
-
- Assert (interpolation_matrix.n() == this->dofs_per_face,
- ExcDimensionMismatch (interpolation_matrix.n(),
- this->dofs_per_face));
- Assert (interpolation_matrix.m() == x_source_fe.dofs_per_face,
- ExcDimensionMismatch (interpolation_matrix.m(),
- x_source_fe.dofs_per_face));
-
- // ok, source is a Q element, so we will be able to do the work
- const FE_Q_Base<POLY,dim,spacedim> &source_fe
- = dynamic_cast<const FE_Q_Base<POLY,dim,spacedim>&>(x_source_fe);
-
- // Make sure that the element for which the DoFs should be constrained is
- // the one with the higher polynomial degree. Actually the procedure will
- // work also if this assertion is not satisfied. But the matrices produced
- // in that case might lead to problems in the hp procedures, which use this
- // method.
- Assert (this->dofs_per_face <= source_fe.dofs_per_face,
- (typename FiniteElement<dim,spacedim>::
- ExcInterpolationNotImplemented ()));
-
- // generate a quadrature with the unit support points. This is later based
- // as a basis for the QProjector, which returns the support points on the
- // face.
- Quadrature<dim-1> quad_face_support (source_fe.get_unit_face_support_points ());
-
- // Rule of thumb for FP accuracy, that can be expected for a given
- // polynomial degree. This value is used to cut off values close to zero.
- const double eps = 2e-13*this->degree*(dim-1);
-
- // compute the interpolation matrix by simply taking the value at the
- // support points.
-//TODO: Verify that all faces are the same with respect to
-// these support points. Furthermore, check if something has to
-// be done for the face orientation flag in 3D.
- const Quadrature<dim> face_quadrature
- = QProjector<dim>::project_to_face (quad_face_support, 0);
- for (unsigned int i=0; i<source_fe.dofs_per_face; ++i)
- {
- const Point<dim> &p = face_quadrature.point (i);
-
- for (unsigned int j=0; j<this->dofs_per_face; ++j)
- {
- double matrix_entry = this->poly_space.compute_value (this->face_to_cell_index(j, 0), p);
-
- // Correct the interpolated value. I.e. if it is close to 1 or 0,
- // make it exactly 1 or 0. Unfortunately, this is required to avoid
- // problems with higher order elements.
- if (std::fabs (matrix_entry - 1.0) < eps)
- matrix_entry = 1.0;
- if (std::fabs (matrix_entry) < eps)
- matrix_entry = 0.0;
-
- interpolation_matrix(i,j) = matrix_entry;
- }
- }
-
- // make sure that the row sum of each of the matrices is 1 at this
- // point. this must be so since the shape functions sum up to 1
- for (unsigned int j=0; j<source_fe.dofs_per_face; ++j)
- {
- double sum = 0.;
-
- for (unsigned int i=0; i<this->dofs_per_face; ++i)
- sum += interpolation_matrix(j,i);
-
- Assert (std::fabs(sum-1) < eps, ExcInternalError());
- }
+ get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int,
+ interpolation_matrix);
}
// these support points. Furthermore, check if something has to
// be done for the face orientation flag in 3D.
const Quadrature<dim> subface_quadrature
- = QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
+ = subface == numbers::invalid_unsigned_int
+ ?
+ QProjector<dim>::project_to_face (quad_face_support, 0)
+ :
+ QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
{
const Point<dim> &p = subface_quadrature.point (i);
}
}
- Assert (best_cell.first.state() == IteratorState::valid,
+ AssertThrow (best_cell.first.state() == IteratorState::valid,
ExcPointNotFound<spacedim>(p));
return best_cell;
}
}
- Assert (best_cell.first.state() == IteratorState::valid,
+ AssertThrow (best_cell.first.state() == IteratorState::valid,
ExcPointNotFound<spacedim>(p));
return best_cell;
}
-
/*
* Internally used in orthogonal_equality
*
*
* See the comment on the next function as well as the detailed
* documentation of make_periodicity_constraints and
- * collect_periodic_face_pairs for details
+ * collect_periodic_faces for details
*/
template<int dim> struct OrientationLookupTable {};
/*
- * Internally used in collect_periodic_face_pairs
+ * Internally used in collect_periodic_faces
+ */
+ template<typename CellIterator>
+ std::vector<PeriodicFacePair<CellIterator> >
+ match_periodic_face_pairs
+ (std::set<std::pair<CellIterator, unsigned int> >
+ &pairs1,
+ std::set<std::pair<typename identity<CellIterator>::type, unsigned int> >
+ &pairs2,
+ int direction,
+ const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
+ &offset)
+ {
+ static const int space_dim = CellIterator::AccessorType::space_dimension;
+ Assert (0<=direction && direction<space_dim,
+ ExcIndexRange (direction, 0, space_dim));
+
+ Assert (pairs1.size() == pairs2.size(),
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ typename std::vector<PeriodicFacePair<CellIterator> > matched_faces;
+
+ // Match with a complexity of O(n^2). This could be improved...
+ std::bitset<3> orientation;
+ typedef typename std::set
+ <std::pair<CellIterator, unsigned int> >::const_iterator PairIterator;
+ for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
+ {
+ for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
+ {
+ const CellIterator cell1 = it1->first;
+ const CellIterator cell2 = it2->first;
+ const unsigned int face_idx1 = it1->second;
+ const unsigned int face_idx2 = it2->second;
+ if (GridTools::orthogonal_equality(orientation,
+ cell1->face(face_idx1),
+ cell2->face(face_idx2),
+ direction, offset))
+ {
+ // We have a match, so insert the matching pairs and
+ // remove the matched cell in pairs2 to speed up the
+ // matching:
+ const PeriodicFacePair<CellIterator> matched_face
+ = {{cell1, cell2},{face_idx1, face_idx2}, orientation};
+ matched_faces.push_back(matched_face);
+ pairs2.erase(it2);
+ break;
+ }
+ }
+ }
+
+ AssertThrow (matched_faces.size() == pairs1.size() && pairs2.size() == 0,
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ return matched_faces;
+ }
+
+ /* Deprecated version of the function above with different return value.
+ * It is used the deprecated collect_periodic_face_pairs.
*/
template<typename FaceIterator>
std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
- match_periodic_face_pairs (std::set<FaceIterator> &faces1, /* not const! */
- std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
- int direction,
- const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ match_periodic_face_pairs
+ (std::set<FaceIterator> &faces1, /* not const! */
+ std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
+ int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
+ DEAL_II_DEPRECATED;
+
+ template<typename FaceIterator>
+ std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
+ match_periodic_face_pairs
+ (std::set<FaceIterator> &faces1, /* not const! */
+ std::set<typename identity<FaceIterator>::type> &faces2, /* not const! */
+ int direction,
+ const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset)
{
static const int space_dim = FaceIterator::AccessorType::space_dimension;
Assert (0<=direction && direction<space_dim,
}
+ template<typename DH>
+ std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ collect_periodic_faces
+ (const DH &dof_handler,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const unsigned int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset)
+ {
+ static const unsigned int dim = DH::dimension;
+ static const unsigned int space_dim = DH::space_dimension;
+ Assert (0<=direction && direction<space_dim,
+ ExcIndexRange (direction, 0, space_dim));
+
+ // Loop over all cells on the highest level and collect all boundary
+ // faces belonging to b_id1 and b_id2:
+
+ std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs1;
+ std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs2;
+
+ for (typename DH::cell_iterator cell = dof_handler.begin(0);
+ cell != dof_handler.end(0); ++cell)
+ {
+ for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
+ {
+ const typename DH::face_iterator face = cell->face(i);
+ if (face->at_boundary() && face->boundary_indicator() == b_id1)
+ {
+ const std::pair<typename DH::cell_iterator, unsigned int> pair1
+ = std::make_pair(cell, i);
+ pairs1.insert(pair1);
+ }
+
+ if (face->at_boundary() && face->boundary_indicator() == b_id2)
+ {
+ const std::pair<typename DH::cell_iterator, unsigned int> pair2
+ = std::make_pair(cell, i);
+ pairs2.insert(pair2);
+ }
+ }
+ }
+
+ Assert (pairs1.size() == pairs2.size(),
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ // and call match_periodic_face_pairs that does the actual matching:
+ return match_periodic_face_pairs(pairs1, pairs2, direction, offset);
+ }
+
+ template<typename DH>
+ std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ collect_periodic_faces (const DH &dof_handler,
+ const types::boundary_id b_id,
+ const unsigned int direction,
+ const dealii::Tensor<1,DH::space_dimension> &offset)
+ {
+ static const unsigned int dim = DH::dimension;
+ static const unsigned int space_dim = DH::space_dimension;
+ Assert (0<=direction && direction<space_dim,
+ ExcIndexRange (direction, 0, space_dim));
+
+ Assert(dim == space_dim,
+ ExcNotImplemented());
+
+ // Loop over all cells on the highest level and collect all boundary
+ // faces 2*direction and 2*direction*1:
+
+ std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs1;
+ std::set<std::pair<typename DH::cell_iterator, unsigned int> > pairs2;
+
+ for (typename DH::cell_iterator cell = dof_handler.begin(0);
+ cell != dof_handler.end(0); ++cell)
+ {
+ const typename DH::face_iterator face_1 = cell->face(2*direction);
+ const typename DH::face_iterator face_2 = cell->face(2*direction+1);
+
+ if (face_1->at_boundary() && face_1->boundary_indicator() == b_id)
+ {
+ const std::pair<typename DH::cell_iterator, unsigned int> pair1
+ = std::make_pair(cell, 2*direction);
+ pairs1.insert(pair1);
+ }
+
+ if (face_2->at_boundary() && face_2->boundary_indicator() == b_id)
+ {
+ const std::pair<typename DH::cell_iterator, unsigned int> pair2
+ = std::make_pair(cell, 2*direction+1);
+ pairs2.insert(pair2);
+ }
+ }
+
+ Assert (pairs1.size() == pairs2.size(),
+ ExcMessage ("Unmatched faces on periodic boundaries"));
+
+ // and call match_periodic_face_pairs that does the actual matching:
+
+ typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ FaceVector;
+
+ FaceVector matching = match_periodic_face_pairs(pairs1, pairs2,
+ direction, offset);
+
+ for (typename FaceVector::iterator pairing = matching.begin();
+ pairing != matching.end(); ++pairing)
+ {
+ Assert(pairing->orientation == 1,
+ ExcMessage("Found a face match with non standard orientation. "
+ "This function is only suitable for meshes with cells "
+ "in default orientation"));
+ }
+
+ return matching;
+ }
template<typename FaceIterator>
std::map<FaceIterator, std::pair<FaceIterator, std::bitset<3> > >
return match_periodic_face_pairs(faces1, faces2, direction, offset);
}
-
-
template<typename DH>
std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > >
collect_periodic_face_pairs (const DH &dof_handler,
int direction,
const dealii::Tensor<1,DH::space_dimension> &offset)
{
- static const int dim = DH::dimension;
- static const int space_dim = DH::space_dimension;
- Assert (0<=direction && direction<space_dim,
- ExcIndexRange (direction, 0, space_dim));
-
- // Loop over all cells on the highest level and collect all boundary
- // faces belonging to b_id1 and b_id2:
+ typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> > FaceVector;
- std::set<typename DH::face_iterator> faces1;
- std::set<typename DH::face_iterator> faces2;
+ const FaceVector face_vector
+ = collect_periodic_faces (dof_handler, b_id1, b_id2, direction, offset);
- for (typename DH::cell_iterator cell = dof_handler.begin(0);
- cell != dof_handler.end(0); ++cell)
- {
- for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
- {
- const typename DH::face_iterator face = cell->face(i);
-
- if (face->at_boundary() && face->boundary_indicator() == b_id1)
- faces1.insert(face);
-
- if (face->at_boundary() && face->boundary_indicator() == b_id2)
- faces2.insert(face);
- }
- }
-
- Assert (faces1.size() == faces2.size(),
- ExcMessage ("Unmatched faces on periodic boundaries"));
+ std::map<typename DH::face_iterator,
+ std::pair<typename DH::face_iterator, std::bitset<3> > >
+ return_value;
+ for(typename FaceVector::const_iterator it = face_vector.begin();
+ it != face_vector.end(); ++it)
+ {
+ const typename DH::face_iterator face1 = it->cell[0]->face(it->face_idx[0]);
+ const typename DH::face_iterator face2 = it->cell[1]->face(it->face_idx[1]);
+ return_value[face1] = std::make_pair(face2, it->orientation);
+ }
- // and call match_periodic_face_pairs that does the actual matching:
- return match_periodic_face_pairs(faces1, faces2, direction, offset);
+ return return_value;
}
-
-
+
template<typename DH>
std::map<typename DH::face_iterator, typename DH::face_iterator>
collect_periodic_face_pairs (const DH &dof_handler,
int direction,
const dealii::Tensor<1,DH::space_dimension> &offset)
{
- static const int dim = DH::dimension;
- static const int space_dim = DH::space_dimension;
- Assert (0<=direction && direction<space_dim,
- ExcIndexRange (direction, 0, space_dim));
-
- Assert(dim == space_dim,
- ExcNotImplemented());
-
- // Loop over all cells on the highest level and collect all boundary
- // faces 2*direction and 2*direction*1:
-
- std::set<typename DH::face_iterator> faces1;
- std::set<typename DH::face_iterator> faces2;
-
- for (typename DH::cell_iterator cell = dof_handler.begin(0);
- cell != dof_handler.end(0); ++cell)
- {
- const typename DH::face_iterator face_1 = cell->face(2*direction);
- const typename DH::face_iterator face_2 = cell->face(2*direction+1);
-
- if (face_1->at_boundary() && face_1->boundary_indicator() == b_id)
- faces1.insert(face_1);
-
- if (face_2->at_boundary() && face_2->boundary_indicator() == b_id)
- faces2.insert(face_2);
- }
-
- Assert (faces1.size() == faces2.size(),
- ExcMessage ("Unmatched faces on periodic boundaries"));
-
- // and call match_periodic_face_pairs that does the actual matching:
-
- typedef std::map<typename DH::face_iterator, std::pair<typename DH::face_iterator, std::bitset<3> > > FaceMap;
- FaceMap matching = match_periodic_face_pairs(faces1, faces2, direction, offset);
-
- std::map<typename DH::face_iterator, typename DH::face_iterator>
- return_value;
-
- for (typename FaceMap::iterator pairing = matching.begin();
- pairing != matching.end(); ++pairing)
- {
- Assert(pairing->second.second == 1,
- ExcMessage("Found a face match with non standard orientation. "
- "This function is only suitable for meshes with cells "
- "in default orientation"));
-
- return_value[pairing->first] = pairing->second.first;
- }
-
+ typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> > FaceVector;
+
+ const FaceVector face_vector
+ = collect_periodic_faces (dof_handler, b_id, direction, offset);
+
+ std::map<typename DH::face_iterator, typename DH::face_iterator> return_value;
+ for(typename FaceVector::const_iterator it = face_vector.begin();
+ it != face_vector.end(); ++it)
+ {
+ const typename DH::face_iterator face1 = it->cell[0]->face(it->face_idx[0]);
+ const typename DH::face_iterator face2 = it->cell[1]->face(it->face_idx[1]);
+ return_value[face1] = face2;
+ }
+
return return_value;
}
-
+
template<typename DH>
void
identify_periodic_face_pairs
typename DH::cell_iterator, unsigned int> >
&periodicity_vector)
{
- static const unsigned int dim = DH::dimension;
- static const unsigned int spacedim = DH::space_dimension;
-
- Assert (0<=direction && direction<spacedim,
- ExcIndexRange(direction, 0, spacedim));
-
- Assert (b_id1 != b_id2,
- ExcMessage ("The boundary indicators b_id1 and b_id2 must be"
- "different to denote different boundaries."));
-
- std::map<std::pair<typename DH::cell_iterator,
- unsigned int>,
- Point<spacedim> > face_locations;
-
- // Collect faces with boundary_indicator b_id1
- typename DH::cell_iterator cell = dof_handler.begin();
- typename DH::cell_iterator endc = dof_handler.end();
- for (; cell != endc; ++cell)
- for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
- if(cell->face(f)->boundary_indicator() == b_id1)
- {
- Point<spacedim> face_center (cell->face(f)->center());
- face_center(direction)=0.;
- const std::pair<typename DH::cell_iterator, unsigned int>
- cell_face_pair = std::make_pair(cell, f);
- face_locations[cell_face_pair]=face_center;
- }
-
- // Match faces with boundary_indicator b_id2 to the ones in
- //face_locations
- cell = dof_handler.begin();
- for (; cell != endc; ++cell)
- for(unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
- if(cell->face(f)->boundary_indicator() == b_id2)
- {
- typename std::map<std::pair<typename DH::cell_iterator,
- unsigned int>,
- Point<spacedim> >::iterator p
- = face_locations.begin();
-
- Point<spacedim> center2 (cell->face(f)->center());
- center2(direction)=0.;
-
- for (; p != face_locations.end(); ++p)
- {
- if (center2.distance(p->second) < 1e-4*cell->face(f)->diameter())
- {
- const std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
- typename DH::cell_iterator, unsigned int>
- periodic_tuple (p->first.first,
- p->first.second,
- cell,
- f);
-
- periodicity_vector.push_back(periodic_tuple);
-
- face_locations.erase(p);
- break;
- }
- Assert (p != face_locations.end(),
- ExcMessage ("No corresponding face was found!"));
- }
- }
+ typedef std::vector<PeriodicFacePair<typename DH::cell_iterator> >
+ FaceVector;
+ const FaceVector periodic_faces
+ = collect_periodic_faces(dof_handler, b_id1, b_id2, direction,
+ dealii::Tensor<1,DH::space_dimension> ());
+
+ typename FaceVector::const_iterator it, end_faces;
+ it = periodic_faces.begin();
+ end_faces = periodic_faces.end();
+
+ for(; it!=end_faces; ++it)
+ {
- Assert (face_locations.size() == 0,
- ExcMessage ("There are unmatched faces!"));
+ const std_cxx1x::tuple<typename DH::cell_iterator, unsigned int,
+ typename DH::cell_iterator, unsigned int>
+ periodic_tuple (it->cell[0], it->face_idx[0],
+ it->cell[1], it->face_idx[1]);
+
+ periodicity_vector.push_back(periodic_tuple);
+ }
}
#if deal_II_dimension >= 2
+ template
+ std::vector<PeriodicFacePair<X::cell_iterator> >
+ collect_periodic_faces
+ (const X &,
+ const types::boundary_id,
+ const types::boundary_id,
+ unsigned int,
+ const Tensor<1,deal_II_space_dimension> &);
+
+ template<typename DH>
+ std::vector<PeriodicFacePair<X::cell_iterator> >
+ collect_periodic_faces
+ (const X &,
+ const types::boundary_id,
+ unsigned int,
+ const Tensor<1,deal_II_space_dimension> &);
+
template
std::map<X::face_iterator, std::pair<X::face_iterator, std::bitset<3> > >
collect_periodic_face_pairs
const types::boundary_id,
int,
const Tensor<1,deal_II_space_dimension> &);
-
template
std::map<X::active_face_iterator, std::pair<X::active_face_iterator, std::bitset<3> > >
collect_periodic_face_pairs
const types::boundary_id,
int,
const Tensor<1,deal_II_space_dimension> &);
-
+
template
void
identify_periodic_face_pairs
// line_dofs_offsets field
dof_handler.faces->lines.dof_offsets
.resize (dof_handler.tria->n_raw_lines(),
- numbers::invalid_dof_index);
+ numbers::invalid_unsigned_int);
unsigned int line_slots_needed = 0;
for (unsigned int line=0; line<dof_handler.tria->n_raw_lines(); ++line)
projects / dealii-svn.git / commitdiff