From 46fc4988d078873fd2275775323fe2d3228165e8 Mon Sep 17 00:00:00 2001 From: heltai Date: Fri, 20 Sep 2013 21:50:33 +0000 Subject: [PATCH] Merged from trunk. git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@30875 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/doxygen/headers/constraints.h | 10 +- deal.II/doc/news/changes.h | 14 + deal.II/doc/publications/index.html | 65 +- deal.II/examples/step-51/doc/intro.dox | 255 +- deal.II/examples/step-51/doc/results.dox | 256 +- deal.II/examples/step-51/step-51.cc | 2131 +++++++++-------- deal.II/include/deal.II/base/graph_coloring.h | 376 +++ .../include/deal.II/base/polynomial_space.h | 2 +- deal.II/include/deal.II/distributed/tria.h | 29 +- deal.II/include/deal.II/dofs/dof_tools.h | 89 +- deal.II/include/deal.II/fe/fe_face.h | 115 +- deal.II/include/deal.II/grid/grid_tools.h | 307 +-- deal.II/include/deal.II/hp/dof_faces.h | 2 +- deal.II/include/deal.II/hp/dof_level.h | 4 +- deal.II/include/deal.II/lac/schur_matrix.h | 2 +- deal.II/source/base/data_out_base.cc | 24 +- deal.II/source/base/mpi.cc | 48 +- deal.II/source/distributed/tria.cc | 108 +- deal.II/source/dofs/dof_tools_constraints.cc | 112 +- .../source/dofs/dof_tools_constraints.inst.in | 8 + deal.II/source/fe/fe_dgp.cc | 49 +- deal.II/source/fe/fe_face.cc | 306 ++- deal.II/source/fe/fe_face.inst.in | 3 +- deal.II/source/fe/fe_q_base.cc | 83 +- deal.II/source/grid/grid_tools.cc | 394 +-- deal.II/source/grid/grid_tools.inst.in | 20 +- deal.II/source/hp/dof_handler.cc | 2 +- 27 files changed, 3004 insertions(+), 1810 deletions(-) create mode 100644 deal.II/include/deal.II/base/graph_coloring.h diff --git a/deal.II/doc/doxygen/headers/constraints.h b/deal.II/doc/doxygen/headers/constraints.h index 999a97659e..0925516053 100644 --- a/deal.II/doc/doxygen/headers/constraints.h +++ b/deal.II/doc/doxygen/headers/constraints.h @@ -26,7 +26,7 @@ * 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 @@ -113,6 +113,13 @@ * 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 positive 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 negative definite. Otherwise, iterative solvers + * such as CG will break down or be considerably slower as GMRES. + * * @note While these two ways are equivalent, 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 @@ -121,7 +128,6 @@ * 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. * - * *

Condensing matrices and sparsity patterns

* * As mentioned above, the first way of using constraints is to build linear diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 89ca970443..870ca7974c 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -85,6 +85,20 @@ inconvenience this causes.

Specific improvements

    + +
  1. + Changed: find_active_cell_around_point now throws the exception ExcPointNotFound + if the point is outside the Triangulation. This exception can be caught. +
    + (Timo Heister, 2013/09/18) +
    2. + +
  2. + Changed: we now call MPI_Init_thread instead of MPI_Init. +
    + (Timo Heister, 2013/09/17) +
    3. +
  3. 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 diff --git a/deal.II/doc/publications/index.html b/deal.II/doc/publications/index.html index 4ce3eff9f5..6630cc3ca3 100644 --- a/deal.II/doc/publications/index.html +++ b/deal.II/doc/publications/index.html @@ -246,6 +246,15 @@ submitted, 2013.
    4. +
  4. G. Alam +
    + Fast iterative solution of large scale statistical + inverse problems + +
    + Master Thesis, Stockholm University, 2013. +
    5. +
  5. J. C. Araujo-Cabarcas
    Numerical methods for glacial isostatic adjustment models @@ -254,13 +263,13 @@ MSc thesis, Uppsala University, Uppsala, Sweden, 2013.
    6. -
  6. G. Alam +
  7. P. Argul, G. Ruocci, F. Freddi, K. Benzarti
    - Fast iterative solution of large scale statistical - inverse problems + An improved damage modelling to deal with the variability of + fracture mechanisms in FRP reinforced concrete structures
    - Master Thesis, Stockholm University, 2013. + International Journal of Adhesion and Adhesives 45 , 7-20, 2013.
  8. D. Arndt @@ -271,6 +280,22 @@ Msc thesis, University of Göttingen, Germany, 2013.
    9. +
  9. A. Barker, T. Rees, M. Stoll +
    + A fast solver for an H1 regularized + PDE-constrained optimization problem +
    + submitted, 2013. +
    10. + +
  10. S. Bartels, R. H. Nochetto, A. J. Salgado +
    + A total variation diminishing interpolation operator + and applications +
    + submitted, 2013. +
    11. +
  11. P. Benner, M. Stoll
    Optimal control for Allen-Cahn equations enhanced by model predictive control @@ -437,6 +462,15 @@ Computers & Mathematics with Applications (in press), 2013.
    12. +
  12. F. Freddi, G. Royer-Carfagni +
    + Plastic flow as an energy minimization problem. Numerical + Experiments + +
    + Journal of Elasticity, in press, 2013. +
    13. +
  13. A. Geringer, B. Lenhof, S. Diebels
    Macroscopic modeling of foams: an order-parameter @@ -511,6 +545,19 @@ pp. 46-72, 2013.
    14. +
  14. S. Kim +
    + Analysis of the convected Helmholtz equation with a + uniform mean flow in a waveguide with complete radiation + boundary conditions + +
    + Journal of Mathematical Analysis and Applications, + in press, + + 2013. +
    15. +
  15. N. Klein
    Consistent FE-Analysis of elliptic Variational Inequalities @@ -588,6 +635,16 @@ a monolithic ALE framework MSc Thesis, University of Cape Town, 2013.
    16. +
  16. T. Reis, W. Wollner +
    + Iterative solution of operator Lyapunov equations + arising in heat transfer + +
    + 2013 European Control Conference, Zurich Switzerland, + pp. 41-46, 2013. +
    17. +
  17. T. Richter, T. Wick
    Optimal Control and Parameter Estimation for Stationary Fluid-Structure Interaction Problems diff --git a/deal.II/examples/step-51/doc/intro.dox b/deal.II/examples/step-51/doc/intro.dox index 4aedf0fe90..7da03aa8c1 100644 --- a/deal.II/examples/step-51/doc/intro.dox +++ b/deal.II/examples/step-51/doc/intro.dox @@ -17,22 +17,23 @@ is the large number of globally coupled degrees of freedom that one 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 for each of the adjacent elements, 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 FE_DGP_Monomial 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.

    Reducing the size of the linear system

    -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: - Hybridizable discontinuous Galerkin methods for partial differential + Hybridizable discontinuous Galerkin methods for partial differential equations in continuum mechanics, Journal of Computational Physics, 2012, 231:18, 5955-5988. [DOI]). @@ -41,15 +42,15 @@ this goal by formulating the mathematical problem using Dirichlet-to-Neumann 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:
      -
    1. 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 +
    2. 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.
    3. Use the known skeleton values as Dirichlet data for solving local element-level solutions. This is known as the -'local solver', and is an embarrassingly parallel element-by-element +'local solver', and is an embarrassingly parallel element-by-element solution process.
    @@ -61,13 +62,13 @@ first block and the skeleton variables $\Lambda$ as the second block: \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 U 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 @@ -78,39 +79,42 @@ The coupling to other cells is introduced by the matrices matrix A 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
      -
    1. constructing the Schur complement matrix $D-C A^{-1} B$ and right hand side $G - C A^{-1} F$, +
    2. constructing the Schur complement matrix $D-C A^{-1} B$ and right hand side $G - C A^{-1} F$ locally on each cell + and inserting the contribution into the global trace matrix in the usual way,
    3. solving the Schur complement system for $\Lambda$, and -
    4. solving the equation for U using the second equation which uses $\Lambda$. +
    5. solving for U using the second equation, given $\Lambda$.

    Solution quality and rates of convergence

    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})$.

    HDG applied to the convection-diffusion problem

    -The HDG formulation used for this example is taken from +The HDG formulation used for this example is taken from
    N.C. Nguyen, J. Peraire, B. Cockburn: - An implicit high-order hybridizable discontinuous Galerkin method - for linear convection–diffusion equations, - Journal of Computational Physics, 2009, 228:9, 3232-3254. + An implicit high-order hybridizable discontinuous Galerkin method + for linear convection–diffusion equations, + Journal of Computational Physics, 2009, 228:9, 3232-3254. [DOI] -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} @@ -120,106 +124,203 @@ the above equation as the first order system: \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_{\partial \mathcal{T}} + + \left_{\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 +uh is defined as the space of functions that are tensor +product polynomials of degree p on each cell and discontinuous over the +element boundaries $\mathcal Q_{-p}$, i.e., the space described by +FE_DGQ(p). The space for the gradient or flux variable +qh 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:
    1. 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)$. -
    2. The matrix $B$ represents the local-face coupling. These are the terms - with weighting functions $(\mathbf{v}, w)$ multiplying the skeleton variable - $\lambda_h$. -
    3. The matrix $C$ represents the face-local coupling, which involves the +
    4. 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$. +
    5. The matrix $C$ represents the face-local coupling, which involves the weighting function $\mu$ multiplying the local solutions $(\mathbf{q}_h, u_h)$. -
    6. The matrix $D$ is the face-face coupling; - terms involve both $\mu$ and $\lambda_h$. +
    7. The matrix $D$ is the face-face coupling; + terms involve both $\mu$ and $\hat{u}_h$.

    Post-processing and super-convergence

    +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 p. For this to happen, there are two necessary +ingredients: +
      +
    1. The computed solution gradient $\mathbf{q}_h$ converges at optimal rate, + i.e., $\mathcal{O}(h^{p+1})$. +
    2. The average of the scalar part of the solution, uh, + on each cell K super-converges at rate $\mathcal{O}(h^{p+2})$. +
    + +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 +K 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 p+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.

    Problem specific data

    -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 xi for the exponentials +
      +
    • 1D: $\{x_i\}^1 = \{ -\frac{1}{3}, 0, \frac{1}{3} \}$, +
    • 2D: $\{\mathbf{x}_i\}^2 = \{ (-\frac{1}{2},\frac{1}{2}), + (-\frac{1}{2},-\frac{1}{2}), + (\frac{1}{2},-\frac{1}{2}) + \}$, +
    • 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}) + \}$. +
    + +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 \\ @@ -227,25 +328,15 @@ $\Omega := [-1,1]^d$ (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$.

    Implementation

    -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: +
      +
    • WorkStream to parallelize local solvers. Workstream is already used in + step-32, step-44. +
    • Reconstruct the local DG solution from the trace. +
    • Post-processing the solution for superconvergence. +
    • DataOutFaces for direct output of the global skeleton solution. +
    diff --git a/deal.II/examples/step-51/doc/results.dox b/deal.II/examples/step-51/doc/results.dox index 5cd183d641..f41ab66065 100644 --- a/deal.II/examples/step-51/doc/results.dox +++ b/deal.II/examples/step-51/doc/results.dox @@ -1,9 +1,97 @@

    Results

    +

    Program output

    -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 p=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. + + + + + + + + + + +
    + + + +
    + + + +
    + +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. + + + + + + + + + + +
    + + + +
    + + + +
    + +Finally, we look at the solution for p=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. + + + + + + +
    + + + +
    + +

    Convergence tables

    + +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: @@ -52,13 +140,13 @@ global refinement was performed, also the convergence rates. The quadratic 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 p converges at rate p only, as opposed -to p+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 p converges at rate p only, as +opposed to p+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. @@ -103,11 +191,151 @@ cells dofs val L2 grad L2 val L2-post 110592 5419008 3.482e-05 3.94 3.055e-04 3.95 7.374e-07 5.00 @endcode +

    Comparison with continuous finite elements

    + +

    Results for 2D

    + +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. + + + + + + +
    + + + +
    + +The results in the graphs show that the HDG method is slower than continuous +finite elements at p=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 p=3 and p=6, and about the same quality +for p=1. + + + + + + +
    + + + +
    + +Since the HDG method actually produces results converging as +hp+2, we should compare it to a continuous Galerkin +solution with the same asymptotic convergence behavior, i.e., FE_Q with degree +p+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. + + + + + + +
    + + + +
    + +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 p=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. + +

    Resuls for 3D

    + +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 +u for CG and HDG at order p, 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 p+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 p but Legendre polynomials of complete +degree p. 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. + + + + + + + + + + + + + + +
    + + + +
    + + + +
    + + + +
    + +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 p=3 and p=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.

    Possibilities for improvements

    -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: @@ -182,7 +410,9 @@ improvement makes most sense. 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.
  18. Speed up assembly by pre-assembling parts that do not change from one cell to another (those that do neither contain variable coefficients nor diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index d385f69061..c8b7a4ffcd 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -80,10 +80,12 @@ // the simulation. #include + +// We start by putting the class into its own namespace. namespace Step51 { -using namespace dealii; + using namespace dealii; // @sect3{Equation data} // @@ -91,98 +93,98 @@ using namespace dealii; // 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 -class SolutionBase -{ -protected: - static const unsigned int n_source_centers = 3; - static const Point 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 + class SolutionBase + { + protected: + static const unsigned int n_source_centers = 3; + static const Point 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 -const double SolutionBase::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 -class Solution : public Function, - protected SolutionBase -{ -public: - Solution () : Function() {} + template + const double SolutionBase::width = 1./5.; - virtual double value (const Point &p, - const unsigned int component = 0) const; - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component = 0) const; -}; + template + class Solution : public Function, + protected SolutionBase + { + public: + Solution () : Function() {} + virtual double value (const Point &p, + const unsigned int component = 0) const; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + }; -template -double Solution::value (const Point &p, - const unsigned int) const -{ - double return_value = 0; - for (unsigned int i=0; in_source_centers; ++i) - { - const Point 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(std::sqrt(2. * numbers::PI) * this->width); -} + template + double Solution::value (const Point &p, + const unsigned int) const + { + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point 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(std::sqrt(2. * numbers::PI) * this->width); + } -template -Tensor<1,dim> Solution::gradient (const Point &p, - const unsigned int) const -{ - Tensor<1,dim> return_value; - for (unsigned int i=0; in_source_centers; ++i) - { - const Point 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 + Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const + { + Tensor<1,dim> return_value; - return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) * - this->width); -} + for (unsigned int i=0; in_source_centers; ++i) + { + const Point 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(std::sqrt(2 * numbers::PI) * + this->width); + } @@ -190,69 +192,69 @@ Tensor<1,dim> Solution::gradient (const Point &p, // 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 -class SolutionAndGradient : public Function, - protected SolutionBase -{ -public: - SolutionAndGradient () : Function(dim) {} + template + class SolutionAndGradient : public Function, + protected SolutionBase + { + public: + SolutionAndGradient () : Function(dim) {} - virtual void vector_value (const Point &p, - Vector &v) const; -}; + virtual void vector_value (const Point &p, + Vector &v) const; + }; -template -void SolutionAndGradient::vector_value (const Point &p, - Vector &v) const -{ - AssertDimension(v.size(), dim+1); - Solution solution; - Tensor<1,dim> grad = solution.gradient(p); - for (unsigned int d=0; d + void SolutionAndGradient::vector_value (const Point &p, + Vector &v) const + { + AssertDimension(v.size(), dim+1); + Solution solution; + Tensor<1,dim> grad = solution.gradient(p); + for (unsigned int d=0; d -class ConvectionVelocity : public TensorFunction<1,dim> -{ -public: - ConvectionVelocity() : TensorFunction<1,dim>() {} + template + class ConvectionVelocity : public TensorFunction<1,dim> + { + public: + ConvectionVelocity() : TensorFunction<1,dim>() {} - virtual Tensor<1,dim> value (const Point &p) const; -}; + virtual Tensor<1,dim> value (const Point &p) const; + }; -template -Tensor<1,dim> -ConvectionVelocity::value(const Point &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 + Tensor<1,dim> + ConvectionVelocity::value(const Point &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; + } @@ -261,42 +263,42 @@ ConvectionVelocity::value(const Point &p) const // 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 -class RightHandSide : public Function, - protected SolutionBase -{ -public: - RightHandSide () : Function() {} + template + class RightHandSide : public Function, + protected SolutionBase + { + public: + RightHandSide () : Function() {} - virtual double value (const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; -private: - const ConvectionVelocity convection_velocity; -}; + private: + const ConvectionVelocity convection_velocity; + }; -template -double RightHandSide::value (const Point &p, - const unsigned int) const -{ - Tensor<1,dim> convection = convection_velocity.value(p); - double return_value = 0; - for (unsigned int i=0; in_source_centers; ++i) - { - const Point 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 + double RightHandSide::value (const Point &p, + const unsigned int) const + { + Tensor<1,dim> convection = convection_velocity.value(p); + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point 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(std::sqrt(2 * numbers::PI) - * this->width); -} + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) + * this->width); + } // @sect3{The HDG solver class} @@ -311,675 +313,683 @@ double RightHandSide::value (const Point &p, // 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 -class HDG -{ -public: - enum RefinementMode + template + 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::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data); + + void copy_local_to_global(const PerTaskData &data); + + void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + unsigned int &empty_data); + + + Triangulation 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 fe_local; + DoFHandler dof_handler_local; + Vector solution_local; + + // The new finite element type and corresponding DoFHandler are + // used for the global skeleton solution that couples the element-level local + // solutions. + FE_FaceQ fe; + DoFHandler dof_handler; + Vector solution; + Vector 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 fe_u_post; + DoFHandler dof_handler_u_post; + Vector 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 ConstrainMatrix 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 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::active_cell_iterator &cell, - ScratchData &scratch, - PerTaskData &task_data); - - void copy_local_to_global(const PerTaskData &data); - - void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, - PostProcessScratchData &scratch, - unsigned int &empty_data); - - - Triangulation 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 fe_local; - DoFHandler dof_handler_local; - Vector solution_local; - - // The new finite element type and corresponding DoFHandler are - // used for the global skeleton solution that couples the element-level local - // solutions. - FE_FaceQ fe; - DoFHandler dof_handler; - Vector solution; - Vector 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 fe_u_post; - DoFHandler dof_handler_u_post; - Vector 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 ConstrainMatrix 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 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 DoFHandler and -// FiniteElement 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 -HDG::HDG (const unsigned int degree, - const RefinementMode refinement_mode) : - fe_local (FE_DGQ(degree), dim, - FE_DGQ(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 DoFHandler objects. The @p solution and @p system_matrix -// objects go with the global skeleton solution. -template -void -HDG::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::type boundary_functions; - Solution solution_function; - boundary_functions[0] = &solution_function; - VectorTools::project_boundary_values (dof_handler, - boundary_functions, - QGauss(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 DoFHandler and + // FiniteElement 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 + HDG::HDG (const unsigned int degree, + const RefinementMode refinement_mode) : + fe_local (FE_DGQ(degree), dim, + FE_DGQ(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 DoFHandler objects. The @p solution and @p system_matrix + // objects go with the global skeleton solution. + template + void + HDG::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::type boundary_functions; + Solution solution_function; + boundary_functions[0] = &solution_function; + VectorTools::project_boundary_values (dof_handler, + boundary_functions, + QGauss(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 -struct HDG::PerTaskData -{ - FullMatrix cell_matrix; - Vector cell_vector; - std::vector 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 WorkStream. The FEValues, matrix, -// and vector objects should be familiar by now. There are two objects that -// need to be discussed: @p std::vector > -// fe_local_support_on_face and @p std::vector > -// 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 -struct HDG::ScratchData -{ - FEValues fe_values_local; - FEFaceValues fe_face_values_local; - FEFaceValues fe_face_values; - - FullMatrix ll_matrix; - FullMatrix lf_matrix; - FullMatrix fl_matrix; - FullMatrix tmp_matrix; - Vector l_rhs; - Vector tmp_rhs; - - std::vector > q_phi; - std::vector q_phi_div; - std::vector u_phi; - std::vector > u_phi_grad; - std::vector tr_phi; - std::vector trace_values; - - std::vector > fe_local_support_on_face; - std::vector > fe_support_on_face; - - ConvectionVelocity convection_velocity; - RightHandSide right_hand_side; - const Solution exact_solution; - - ScratchData(const FiniteElement &fe, - const FiniteElement &fe_local, - const QGauss &quadrature_formula, - const QGauss &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::faces_per_cell), - fe_support_on_face(GeometryInfo::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 + struct HDG::PerTaskData { - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int i=0; i cell_matrix; + Vector cell_vector; + std::vector 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::faces_per_cell; ++face) - for (unsigned int i=0; iWorkStream. The FEValues, matrix, + // and vector objects should be familiar by now. There are two objects that + // need to be discussed: @p std::vector > + // fe_local_support_on_face and @p std::vector > + // 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 + struct HDG::ScratchData + { + FEValues fe_values_local; + FEFaceValues fe_face_values_local; + FEFaceValues fe_face_values; + + FullMatrix ll_matrix; + FullMatrix lf_matrix; + FullMatrix fl_matrix; + FullMatrix tmp_matrix; + Vector l_rhs; + Vector tmp_rhs; + + std::vector > q_phi; + std::vector q_phi_div; + std::vector u_phi; + std::vector > u_phi_grad; + std::vector tr_phi; + std::vector trace_values; + + std::vector > fe_local_support_on_face; + std::vector > fe_support_on_face; + + ConvectionVelocity convection_velocity; + RightHandSide right_hand_side; + const Solution exact_solution; + + ScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const QGauss &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::faces_per_cell), + fe_support_on_face(GeometryInfo::faces_per_cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int i=0; i::faces_per_cell; ++face) + for (unsigned int i=0; iWorkStream -// when post-processing the local solution $u^*$. It is similar, but much -// simpler, than @p ScratchData. -template -struct HDG::PostProcessScratchData -{ - FEValues fe_values_local; - FEValues fe_values; - - std::vector u_values; - std::vector > u_gradients; - FullMatrix cell_matrix; - - Vector cell_rhs; - Vector cell_sol; - - PostProcessScratchData(const FiniteElement &fe, - const FiniteElement &fe_local, - const QGauss &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 WorkStream + // when post-processing the local solution $u^*$. It is similar, but much + // simpler, than @p ScratchData. + template + struct HDG::PostProcessScratchData + { + FEValues fe_values_local; + FEValues fe_values; + + std::vector u_values; + std::vector > u_gradients; + FullMatrix cell_matrix; + + Vector cell_rhs; + Vector cell_sol; + + PostProcessScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &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 -void HDG::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 Step-32, where + // the quadrature formula and the update flags are set up, and then + // WorkStream 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 + void + HDG::assemble_system (const bool trace_reconstruct) + { + const QGauss quadrature_formula(fe.degree+1); + const QGauss 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 Step-32, where -// the quadrature formula and the update flags are set up, and then -// WorkStream 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 -void -HDG::assemble_system (const bool trace_reconstruct) -{ - const QGauss quadrature_formula(fe.degree+1); - const QGauss 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::assemble_system_one_cell, - &HDG::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::assemble_system_one_cell, + &HDG::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 -void -HDG::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, - ScratchData &scratch, - PerTaskData &task_data) -{ -// Construct iterator for dof_handler_local for FEValues reinit function. - typename DoFHandler::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 + void + HDG::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data) + { + // Construct iterator for dof_handler_local for FEValues reinit function. + typename DoFHandler::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 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::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 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 quadrature_point = - scratch.fe_face_values.quadrature_point(q); - const Point 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::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 quadrature_point = + scratch.fe_face_values.quadrature_point(q); + const Point 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; kface(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_{\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; iface(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_{\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; iget_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 set_dof_values 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 set_dof_values 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 -void HDG::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 + void HDG::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 + void HDG::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); + } @@ -1005,69 +1015,69 @@ void HDG::solve () // 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 -void -HDG::postprocess() -{ + template + void + HDG::postprocess() { - const QGauss 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 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::postprocess_one_cell, + std_cxx1x::ref(*this), + std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), + std_cxx1x::function(), + 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::postprocess_one_cell, - std_cxx1x::ref(*this), - std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), - std_cxx1x::function(), - scratch, - 0U); + Vector difference_per_cell (triangulation.n_active_cells()); + + ComponentSelectFunction value_select (dim, dim+1); + VectorTools::integrate_difference (dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &value_select); + const double L2_error = difference_per_cell.l2_norm(); + + ComponentSelectFunction gradient_select (std::pair(0, dim), + dim+1); + VectorTools::integrate_difference (dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(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(), + difference_per_cell, + QGauss(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 difference_per_cell (triangulation.n_active_cells()); - - ComponentSelectFunction value_select (dim, dim+1); - VectorTools::integrate_difference (dof_handler_local, - solution_local, - SolutionAndGradient(), - difference_per_cell, - QGauss(fe.degree+2), - VectorTools::L2_norm, - &value_select); - const double L2_error = difference_per_cell.l2_norm(); - - ComponentSelectFunction gradient_select (std::pair(0, dim), - dim+1); - VectorTools::integrate_difference (dof_handler_local, - solution_local, - SolutionAndGradient(), - difference_per_cell, - QGauss(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(), - difference_per_cell, - QGauss(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} @@ -1076,7 +1086,7 @@ HDG::postprocess() // 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 @@ -1089,154 +1099,154 @@ HDG::postprocess() // 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 -void -HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, - PostProcessScratchData &scratch, - unsigned int &) -{ - typename DoFHandler::active_cell_iterator - loc_cell (&triangulation, - cell->level(), - cell->index(), - &dof_handler_local); + template + void + HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + unsigned int &) + { + typename DoFHandler::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; idistribute_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 DoFHandler -// objects. The graphical output for the skeleton variable is done through -// use of the DataOutFaces class. -template -void HDG::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 DoFHandler + // objects. The graphical output for the skeleton variable is done through + // use of the DataOutFaces class. + template + void HDG::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 data_out; - -// We first define the names and types of the local solution, -// and add the data to @p data_out. - std::vector names (dim, "gradient"); - names.push_back ("solution"); - std::vector - 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 post_name(1,"u_post"); - std::vector - 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 data_out; + + // We first define the names and types of the local solution, + // and add the data to @p data_out. + std::vector names (dim, "gradient"); + names.push_back ("solution"); + std::vector + 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 post_name(1,"u_post"); + std::vector + 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 DataOutFaces class works analagously to the DataOut // class when we have a DoFHandler 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 data_out_face(false); - std::vector face_name(1,"lambda"); - std::vector - 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 data_out_face(false); + std::vector face_name(1,"u_hat"); + std::vector + 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} @@ -1250,119 +1260,120 @@ void HDG::output_results (const unsigned int cycle) // The adaptive_refinement mode uses the KellyErrorEstimator to // give a decent indication of the non-regular regions in the scalar local // solutions. -template -void HDG::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 + void HDG::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 estimated_error_per_cell (triangulation.n_active_cells()); + case adaptive_refinement: + { + Vector estimated_error_per_cell (triangulation.n_active_cells()); - FEValuesExtractors::Scalar scalar(dim); - typename FunctionMap::type neumann_boundary; - KellyErrorEstimator::estimate (dof_handler_local, - QGauss(3), - neumann_boundary, - solution_local, - estimated_error_per_cell, - fe_local.component_mask(scalar)); + FEValuesExtractors::Scalar scalar(dim); + typename FunctionMap::type neumann_boundary; + KellyErrorEstimator::estimate (dof_handler_local, + QGauss(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::cell_iterator + cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::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 Step-7. + // We loop over 10 cycles, refining the grid on each one. At the end, + // convergence tables are created. + template + void HDG::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::cell_iterator - cell = triangulation.begin (), - endc = triangulation.end(); - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::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 Step-7. -// We loop over 10 cycles, refining the grid on each one. At the end, -// convergence tables are created. -template -void HDG::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) { diff --git a/deal.II/include/deal.II/base/graph_coloring.h b/deal.II/include/deal.II/base/graph_coloring.h new file mode 100644 index 0000000000..3bf6c6abf5 --- /dev/null +++ b/deal.II/include/deal.II/base/graph_coloring.h @@ -0,0 +1,376 @@ + +// --------------------------------------------------------------------- +// $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 +#include +#include +#include + +#include + +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 + std::vector > create_partitioning(Iterator const &begin, + typename identity::type const &end, + std::vector (*get_conflict_indices) (Iterator &)) + { + std::vector > partitioning(1,std::vector (1,begin)); + // Number of iterators. + unsigned int n_iterators(0); + // Create a map from conflict indices to iterators + boost::unordered_map > indices_to_iterators; + for (Iterator it=begin; it!=end; ++it) + { + std::vector conflict_indices = (*get_conflict_indices)(it); + const unsigned int n_conflict_indices(conflict_indices.size()); + for (unsigned int i=0; i used_it; + used_it.insert(begin); + while (used_it.size()!=n_iterators) + { + typename std::vector::iterator vector_it(partitioning.back().begin()); + typename std::vector::iterator vector_end(partitioning.back().end()); + std::vector new_zone; + for (; vector_it!=vector_end; ++vector_it) + { + std::vector conflict_indices = (*get_conflict_indices)(*vector_it); + const unsigned int n_conflict_indices(conflict_indices.size()); + for (unsigned int i=0; i iterator_vector(indices_to_iterators[conflict_indices[i]]); + for (unsigned int j=0; j (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 + std::vector > make_dsatur_coloring(std::vector &partition, + std::vector (*get_conflict_indices)(Iterator &)) + { + std::vector > partition_coloring; + // Number of zones composing the partitioning. + const unsigned int partition_size(partition.size()); + std::vector sorted_vertices(partition_size); + std::vector degrees(partition_size); + std::vector > conflict_indices(partition_size); + std::vector > 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 conflict_indices_intersection; + std::vector::iterator intersection_it; + for (unsigned int i=0; i::iterator degrees_it; + for (unsigned int i=0; i > colors_used; + for (unsigned int i=0; i (1, + partition[current_vertex])); + boost::unordered_set 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 + std::vector > + gather_colors(std::vector > > const &partition_coloring) + { + std::vector > coloring; + + // Count the number of iterators in each color. + const unsigned int partition_size(partition_coloring.size()); + std::vector > colors_counter(partition_size); + for (unsigned int i=0; i used_k; + for (unsigned int j=0; j::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 used_k; + for (unsigned int j=0; j::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 + std::vector > + make_graph_coloring(Iterator const &begin,typename identity::type const &end, + std::vector (*get_conflict_indices)(Iterator &)) + { + // Create the partitioning. + std::vector > partitioning = create_partitioning(begin,end, + get_conflict_indices); + + // Color the iterators within each partition. + const unsigned int partitioning_size(partitioning.size()); + std::vector > > partition_coloring( + partitioning_size); + for (unsigned int i=0; i > 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 --------------------------- diff --git a/deal.II/include/deal.II/base/polynomial_space.h b/deal.II/include/deal.II/base/polynomial_space.h index 7865f3f581..56a9f26e9e 100644 --- a/deal.II/include/deal.II/base/polynomial_space.h +++ b/deal.II/include/deal.II/base/polynomial_space.h @@ -236,7 +236,7 @@ protected: * pi(x)pj(y)pk(z). */ void compute_index (const unsigned int n, - unsigned int (&index)[dim]) const; + unsigned int (&index)[dim>0?dim:1]) const; private: /** diff --git a/deal.II/include/deal.II/distributed/tria.h b/deal.II/include/deal.II/distributed/tria.h index 95c34afc1d..5d647258f1 100644 --- a/deal.II/include/deal.II/distributed/tria.h +++ b/deal.II/include/deal.II/distributed/tria.h @@ -158,7 +158,8 @@ namespace internal } } - +//forward declaration of the data type for periodic face pairs +namespace GridTools {template struct PeriodicFacePair;} namespace parallel { @@ -700,15 +701,36 @@ 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 >&); + + /** + * Same as the function above, but takes a different argument. + * * The entries in the std::vector should have the form * std_cxx1x::tuple. * * 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: @@ -720,6 +742,7 @@ namespace parallel cell_iterator, unsigned int> >&); + private: /** * MPI communicator to be diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h index e4d237faf3..9b39197104 100644 --- a/deal.II/include/deal.II/dofs/dof_tools.h +++ b/deal.II/include/deal.II/dofs/dof_tools.h @@ -53,6 +53,7 @@ class ConstraintMatrix; template class InterGridMap; template class Mapping; +namespace GridTools {template 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 @@ -1087,13 +1088,15 @@ namespace DoFTools */ template void - make_periodicity_constraints (const FaceIterator &face_1, - const typename identity::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::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); + /** @@ -1143,12 +1146,13 @@ namespace DoFTools */ template 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()); /** @@ -1167,17 +1171,18 @@ namespace DoFTools * * @see @ref GlossBoundaryIndicator "Glossary entry on boundary indicators" * - * @author Matthias Maier, 2012 + * @author Daniel Arndt, 2013, Matthias Maier, 2012 */ template 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()); /** @@ -1204,11 +1209,12 @@ namespace DoFTools */ template 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()); /** @@ -1233,12 +1239,31 @@ namespace DoFTools */ template 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 + void + make_periodicity_constraints + (const std::vector > + &periodic_faces, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask = ComponentMask()); //@} diff --git a/deal.II/include/deal.II/fe/fe_face.h b/deal.II/include/deal.II/fe/fe_face.h index 3047f7a06a..64bb4e4f2b 100644 --- a/deal.II/include/deal.II/fe/fe_face.h +++ b/deal.II/include/deal.II/fe/fe_face.h @@ -19,6 +19,7 @@ #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -32,8 +33,9 @@ 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 @@ -105,10 +107,112 @@ public: /** * 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 &fe_other) const; + +private: + /** + * Return vector with dofs per vertex, line, quad, hex. + */ + static std::vector 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 +class FE_FaceP : public FE_PolyFace, dim, spacedim> +{ +public: + /** + * Constructor for complete basis of polynomials of degree p. The + * shape functions created using this constructor correspond to Legendre + * polynomials in each coordinate direction. + */ + FE_FaceP(unsigned int p); + + /** + * Clone method. + */ + virtual FiniteElement *clone() const; + + /** + * Return a string that uniquely identifies a finite element. This class + * returns FE_FaceP(degree) , with dim and + * degree 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 + * source.dofs_per_face times this->dofs_per_face. This + * element only provides interpolation matrices for elements of the same + * type and FE_Nothing. For all other elements, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. + */ + virtual void + get_face_interpolation_matrix (const FiniteElement &source, + FullMatrix &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 + * source.dofs_per_face times this->dofs_per_face. This + * element only provides interpolation matrices for elements of the same + * type and FE_Nothing. For all other elements, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. + */ + virtual void + get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &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; @@ -131,6 +235,7 @@ private: static std::vector get_dpo_vector (const unsigned int deg); }; + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/deal.II/include/deal.II/grid/grid_tools.h b/deal.II/include/deal.II/grid/grid_tools.h index 83e40cf654..29900ec738 100644 --- a/deal.II/include/deal.II/grid/grid_tools.h +++ b/deal.II/include/deal.II/grid/grid_tools.h @@ -344,6 +344,11 @@ namespace GridTools * 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. @@ -355,32 +360,23 @@ namespace GridTools /** - * 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 @@ -396,58 +392,39 @@ namespace GridTools const Point &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 @@ -464,14 +441,16 @@ namespace GridTools const Point &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 @@ -990,7 +969,18 @@ namespace GridTools const std::set &boundary_ids = std::set()); - + /** + * 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 + struct PeriodicFacePair + { + CellIterator cell[2]; + unsigned int face_idx[2]; + std::bitset<3> orientation; + }; /** * An orthogonal equality test for faces. @@ -1076,8 +1066,8 @@ namespace GridTools /** - * 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 @@ -1087,45 +1077,33 @@ namespace GridTools * 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 - std::map > > - collect_periodic_face_pairs (const FaceIterator &begin, - const typename identity::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 - std::map > > - 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 > + 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 @@ -1145,40 +1123,101 @@ namespace GridTools * meshes with cells not in @ref GlossFaceOrientation * "standard orientation". * + * @author Daniel Arndt, 2013 + */ + template + std::vector > + 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 + std::map > > + 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 std::map - 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 + std::map > > + collect_periodic_face_pairs + (const FaceIterator &begin, + const typename identity::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 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 > - &periodicity_vector); + (const DH &dof_handler, + const types::boundary_id b_id1, + const types::boundary_id b_id2, + const unsigned int direction, + std::vector > + &periodicity_vector) DEAL_II_DEPRECATED; /** diff --git a/deal.II/include/deal.II/hp/dof_faces.h b/deal.II/include/deal.II/hp/dof_faces.h index 0f3f2cbb0e..a5bce7eedc 100644 --- a/deal.II/include/deal.II/hp/dof_faces.h +++ b/deal.II/include/deal.II/hp/dof_faces.h @@ -495,7 +495,7 @@ namespace internal // 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 const MDt &Dt, const MC &C, VectorMemory > &mem, - const std::vector &signature) + const std::vector &signature) : Ainv(&Ainv), B(&B), Dt(&Dt), C(&C), mem(mem), signature(signature), diff --git a/deal.II/source/base/data_out_base.cc b/deal.II/source/base/data_out_base.cc index b285026281..e1238da38b 100644 --- a/deal.II/source/base/data_out_base.cc +++ b/deal.II/source/base/data_out_base.cc @@ -1050,7 +1050,7 @@ namespace */ std::vector vertices; std::vector cells; - unsigned int cell_offset; + unsigned int node_offset; }; @@ -1556,7 +1556,7 @@ namespace 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 @@ -1577,18 +1577,18 @@ namespace unsigned int d3) { unsigned int base_entry = index * GeometryInfo::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; } } } @@ -6768,7 +6768,7 @@ XDMFEntry DataOutBase::create_xdmf_entry (const std::vector 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; @@ -7124,7 +7124,7 @@ void DataOutBase::write_hdf5_parallel (const std::vector > & 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; diff --git a/deal.II/source/base/mpi.cc b/deal.II/source/base/mpi.cc index 663d647dad..4b93f9be56 100644 --- a/deal.II/source/base/mpi.cc +++ b/deal.II/source/base/mpi.cc @@ -342,16 +342,9 @@ namespace Utilities 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; @@ -359,16 +352,36 @@ namespace Utilities 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 @@ -426,11 +439,11 @@ namespace Utilities #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); @@ -451,7 +464,6 @@ namespace Utilities AssertThrow (mpi_err == 0, ExcMessage ("An error occurred while calling MPI_Finalize()")); -#endif #endif } diff --git a/deal.II/source/distributed/tria.cc b/deal.II/source/distributed/tria.cc index ce4cbc6e97..60f6249249 100644 --- a/deal.II/source/distributed/tria.cc +++ b/deal.II/source/distributed/tria.cc @@ -23,8 +23,8 @@ #include #include #include -#include #include +#include #ifdef DEAL_II_WITH_P4EST # include @@ -3436,39 +3436,31 @@ namespace parallel return mpi_communicator; } - - template + template void Triangulation::add_periodicity - (const std::vector >& - periodicity_vector) + (const std::vector >& + 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 - >::const_iterator periodic_tuple; - periodic_tuple = periodicity_vector.begin(); - - typename std::vector - >::const_iterator periodic_end; + + typedef std::vector > + FaceVector; + typename FaceVector::const_iterator it, periodic_end; + it = periodicity_vector.begin(); periodic_end = periodicity_vector.end(); - - for (; periodic_tuple(*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(), @@ -3477,6 +3469,12 @@ namespace parallel = 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:: connectivity_join_faces (connectivity, tree_left, @@ -3484,6 +3482,7 @@ namespace parallel 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, @@ -3495,25 +3494,60 @@ namespace parallel * date. */ } - - + + Assert(dealii::internal::p4est::functions::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::destroy (parallel_forest); parallel_forest = dealii::internal::p4est::functions:: 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 + void + Triangulation::add_periodicity + (const std::vector >& + periodicity_vector) + { +#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1) + typedef std::vector > + FaceVector; + + typename FaceVector::const_iterator it, end_periodic; + it = periodicity_vector.begin(); + end_periodic = periodicity_vector.end(); + + std::vector > 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 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 diff --git a/deal.II/source/dofs/dof_tools_constraints.cc b/deal.II/source/dofs/dof_tools_constraints.cc index 402557bee3..4909d3a7b5 100644 --- a/deal.II/source/dofs/dof_tools_constraints.cc +++ b/deal.II/source/dofs/dof_tools_constraints.cc @@ -1970,7 +1970,7 @@ namespace DoFTools } - + template void make_periodicity_constraints (const DH &dof_handler, @@ -2010,43 +2010,16 @@ namespace DoFTools 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 > > FaceMap; + typedef std::vector > 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 + (matched_faces, constraint_matrix, component_mask); } @@ -2088,40 +2061,57 @@ namespace DoFTools Assert(dim == space_dim, ExcNotImplemented()); - typedef typename DH::face_iterator FaceIterator; - typedef std::map FaceMap; + typedef std::vector > 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 + (matched_faces, constraint_matrix, component_mask); + } + + + + template + void + make_periodicity_constraints + (const std::vector > + &periodic_faces, + dealii::ConstraintMatrix &constraint_matrix, + const ComponentMask &component_mask) + { + typedef std::vector > + 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]); + } } diff --git a/deal.II/source/dofs/dof_tools_constraints.inst.in b/deal.II/source/dofs/dof_tools_constraints.inst.in index 246f53aa02..d61208de26 100644 --- a/deal.II/source/dofs/dof_tools_constraints.inst.in +++ b/deal.II/source/dofs/dof_tools_constraints.inst.in @@ -68,6 +68,14 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS) dealii::Tensor<1,DH::space_dimension> &, dealii::ConstraintMatrix &, const ComponentMask &); + + template + void + DoFTools::make_periodicity_constraints + (const std::vector > &, + dealii::ConstraintMatrix &, + const ComponentMask &); + #endif } diff --git a/deal.II/source/fe/fe_dgp.cc b/deal.II/source/fe/fe_dgp.cc index bb3f9531af..1c4d07bc3e 100644 --- a/deal.II/source/fe/fe_dgp.cc +++ b/deal.II/source/fe/fe_dgp.cc @@ -32,9 +32,8 @@ FE_DGP::FE_DGP (const unsigned int degree) std::vector(FiniteElementData( get_dpo_vector(degree), 1, degree).dofs_per_cell, std::vector(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) @@ -50,12 +49,9 @@ template std::string FE_DGP::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 << ")"; @@ -101,12 +97,10 @@ FE_DGP:: get_face_interpolation_matrix (const FiniteElement &x_source_fe, FullMatrix &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 FE; typedef FE_DGP FEDGP; AssertThrow ((x_source_fe.get_name().find ("FE_DGP<") == 0) @@ -132,12 +126,10 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe const unsigned int , FullMatrix &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 FE; typedef FE_DGP FEDGP; AssertThrow ((x_source_fe.get_name().find ("FE_DGP<") == 0) @@ -170,8 +162,7 @@ std::vector > FE_DGP:: hp_vertex_dof_identities (const FiniteElement &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*>(&fe_other) != 0) return std::vector > (); @@ -189,8 +180,7 @@ std::vector > FE_DGP:: hp_line_dof_identities (const FiniteElement &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*>(&fe_other) != 0) return std::vector > (); @@ -208,8 +198,7 @@ std::vector > FE_DGP:: hp_quad_dof_identities (const FiniteElement &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*>(&fe_other) != 0) return std::vector > (); @@ -226,8 +215,7 @@ template FiniteElementDomination::Domination FE_DGP::compare_for_face_domination (const FiniteElement &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*>(&fe_other) != 0) return FiniteElementDomination::no_requirements; @@ -243,8 +231,7 @@ bool FE_DGP::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; } diff --git a/deal.II/source/fe/fe_face.cc b/deal.II/source/fe/fe_face.cc index bbf8124ada..641e20a684 100644 --- a/deal.II/source/fe/fe_face.cc +++ b/deal.II/source/fe/fe_face.cc @@ -18,7 +18,9 @@ #include #include #include +#include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -73,7 +75,8 @@ FE_FaceQ::FE_FaceQ (const unsigned int degree) 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::faces_per_cell* this->unit_face_support_points.size()); const unsigned int n_face_dofs = this->unit_face_support_points.size(); @@ -94,6 +97,7 @@ FE_FaceQ::FE_FaceQ (const unsigned int degree) } + template FiniteElement * FE_FaceQ::clone() const @@ -102,17 +106,14 @@ FE_FaceQ::clone() const } + template std::string FE_FaceQ::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 << ")"; @@ -124,76 +125,11 @@ FE_FaceQ::get_name () const template void FE_FaceQ:: -get_face_interpolation_matrix (const FiniteElement &x_source_fe, +get_face_interpolation_matrix (const FiniteElement &source_fe, FullMatrix &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 *>(&x_source_fe) != 0), - (typename FiniteElement:: - 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 &source_fe - = dynamic_cast&>(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:: - ExcInterpolationNotImplemented ())); - - // generate a quadrature with the unit face support points. - const Quadrature 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 &p = face_quadrature.point (i); - - for (unsigned int j=0; jdofs_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; jdofs_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); } @@ -241,6 +177,10 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe for (unsigned int i=0; idofs_per_face; ++i) { const Point p = + subface == numbers::invalid_unsigned_int + ? + face_quadrature.point(i) + : GeometryInfo::child_to_cell_coordinates (face_quadrature.point(i), subface); @@ -342,6 +282,220 @@ compare_for_face_domination (const FiniteElement &fe_other) const return FiniteElementDomination::neither_element_dominates; } + + +// --------------------------------------- FE_FaceP -------------------------- + +template +FE_FaceP::FE_FaceP (const unsigned int degree) + : + FE_PolyFace, dim, spacedim> + (PolynomialSpace(Polynomials::Legendre::generate_complete_basis(degree)), + FiniteElementData(get_dpo_vector(degree), 1, degree, FiniteElementData::L2), + std::vector(1,true)) +{} + + +template +FiniteElement * +FE_FaceP::clone() const +{ + return new FE_FaceP(this->degree); +} + + +template +std::string +FE_FaceP::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 +bool +FE_FaceP::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 +std::vector +FE_FaceP::get_dpo_vector (const unsigned int deg) +{ + std::vector dpo(dim+1, 0U); + dpo[dim-1] = deg+1; + for (unsigned int i=1; i +bool +FE_FaceP::hp_constraints_are_implemented () const +{ + return true; +} + + + +template +FiniteElementDomination::Domination +FE_FaceP:: +compare_for_face_domination (const FiniteElement &fe_other) const +{ + if (const FE_FaceP *fe_q_other + = dynamic_cast*>(&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*>(&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 +void +FE_FaceP:: +get_face_interpolation_matrix (const FiniteElement &source_fe, + FullMatrix &interpolation_matrix) const +{ + get_subface_interpolation_matrix (source_fe, numbers::invalid_unsigned_int, + interpolation_matrix); +} + + + +template +void +FE_FaceP:: +get_subface_interpolation_matrix (const FiniteElement &x_source_fe, + const unsigned int subface, + FullMatrix &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 *source_fe + = dynamic_cast *>(&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:: + 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 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 mass (face_quadrature.size(), source_fe->dofs_per_face); + + for (unsigned int k = 0; k < face_quadrature.size(); ++k) + { + const Point p = + subface == numbers::invalid_unsigned_int ? + face_quadrature.point(k) : + GeometryInfo::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 H(mass); + Vector v_in(face_quadrature.size()); + Vector 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; idofs_per_face; ++i) + { + for (unsigned int k = 0; k < face_quadrature.size(); ++k) + { + const Point p = numbers::invalid_unsigned_int ? + face_quadrature.point(k) : + GeometryInfo::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 *>(&x_source_fe) != 0) + { + // nothing to do here, the FE_Nothing has no degrees of freedom anyway + } + else + AssertThrow (false,(typename FiniteElement:: + ExcInterpolationNotImplemented())); +} + + // explicit instantiations #include "fe_face.inst" diff --git a/deal.II/source/fe/fe_face.inst.in b/deal.II/source/fe/fe_face.inst.in index 8d99474428..513eb0030b 100644 --- a/deal.II/source/fe/fe_face.inst.in +++ b/deal.II/source/fe/fe_face.inst.in @@ -20,9 +20,10 @@ for (deal_II_dimension : DIMENSIONS) { #if deal_II_dimension > 1 template class FE_PolyFace >; - //template class FE_PolyFace, deal_II_dimension>; + template class FE_PolyFace, deal_II_dimension>; //template class FE_PolyFace, deal_II_dimension>; template class FE_FaceQ; + template class FE_FaceP; #endif } diff --git a/deal.II/source/fe/fe_q_base.cc b/deal.II/source/fe/fe_q_base.cc index 6a73cb86ff..40096c09cf 100644 --- a/deal.II/source/fe/fe_q_base.cc +++ b/deal.II/source/fe/fe_q_base.cc @@ -550,83 +550,12 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, template void FE_Q_Base:: -get_face_interpolation_matrix (const FiniteElement &x_source_fe, +get_face_interpolation_matrix (const FiniteElement &source_fe, FullMatrix &interpolation_matrix) const { Assert (dim > 1, ExcImpossibleInDim(1)); - - // this is only implemented, if the source FE is also a Q element - AssertThrow ((dynamic_cast *>(&x_source_fe) != 0), - (typename FiniteElement:: - 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 &source_fe - = dynamic_cast&>(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:: - 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 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 face_quadrature - = QProjector::project_to_face (quad_face_support, 0); - for (unsigned int i=0; i &p = face_quadrature.point (i); - - for (unsigned int j=0; jdofs_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; jdofs_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); } @@ -676,7 +605,11 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe // these support points. Furthermore, check if something has to // be done for the face orientation flag in 3D. const Quadrature subface_quadrature - = QProjector::project_to_subface (quad_face_support, 0, subface); + = subface == numbers::invalid_unsigned_int + ? + QProjector::project_to_face (quad_face_support, 0) + : + QProjector::project_to_subface (quad_face_support, 0, subface); for (unsigned int i=0; idofs_per_face; ++i) { const Point &p = subface_quadrature.point (i); diff --git a/deal.II/source/grid/grid_tools.cc b/deal.II/source/grid/grid_tools.cc index eeeb0dc6ad..4f6981479d 100644 --- a/deal.II/source/grid/grid_tools.cc +++ b/deal.II/source/grid/grid_tools.cc @@ -965,7 +965,7 @@ next_cell: } } - Assert (best_cell.first.state() == IteratorState::valid, + AssertThrow (best_cell.first.state() == IteratorState::valid, ExcPointNotFound(p)); return best_cell; @@ -1088,7 +1088,7 @@ next_cell: } } - Assert (best_cell.first.state() == IteratorState::valid, + AssertThrow (best_cell.first.state() == IteratorState::valid, ExcPointNotFound(p)); return best_cell; @@ -2217,7 +2217,6 @@ next_cell: } - /* * Internally used in orthogonal_equality * @@ -2226,7 +2225,7 @@ next_cell: * * 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 struct OrientationLookupTable {}; @@ -2355,14 +2354,82 @@ next_cell: /* - * Internally used in collect_periodic_face_pairs + * Internally used in collect_periodic_faces + */ + template + std::vector > + match_periodic_face_pairs + (std::set > + &pairs1, + std::set::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 > matched_faces; + + // Match with a complexity of O(n^2). This could be improved... + std::bitset<3> orientation; + typedef typename std::set + >::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 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 std::map > > - match_periodic_face_pairs (std::set &faces1, /* not const! */ - std::set::type> &faces2, /* not const! */ - int direction, - const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + match_periodic_face_pairs + (std::set &faces1, /* not const! */ + std::set::type> &faces2, /* not const! */ + int direction, + const dealii::Tensor<1,FaceIterator::AccessorType::space_dimension> &offset) + DEAL_II_DEPRECATED; + + template + std::map > > + match_periodic_face_pairs + (std::set &faces1, /* not const! */ + std::set::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 + std::vector > + 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 > pairs1; + std::set > pairs2; + + for (typename DH::cell_iterator cell = dof_handler.begin(0); + cell != dof_handler.end(0); ++cell) + { + for (unsigned int i = 0; i < GeometryInfo::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 pair1 + = std::make_pair(cell, i); + pairs1.insert(pair1); + } + + if (face->at_boundary() && face->boundary_indicator() == b_id2) + { + const std::pair 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 + std::vector > + 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 > pairs1; + std::set > 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 pair1 + = std::make_pair(cell, 2*direction); + pairs1.insert(pair1); + } + + if (face_2->at_boundary() && face_2->boundary_indicator() == b_id) + { + const std::pair 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 > + 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 std::map > > @@ -2436,8 +2616,6 @@ next_cell: return match_periodic_face_pairs(faces1, faces2, direction, offset); } - - template std::map > > collect_periodic_face_pairs (const DH &dof_handler, @@ -2446,41 +2624,26 @@ next_cell: 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 > FaceVector; - std::set faces1; - std::set 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::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 > > + 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 std::map collect_periodic_face_pairs (const DH &dof_handler, @@ -2488,59 +2651,24 @@ next_cell: 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 faces1; - std::set 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 > > FaceMap; - FaceMap matching = match_periodic_face_pairs(faces1, faces2, direction, offset); - - std::map - 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 > FaceVector; + + const FaceVector face_vector + = collect_periodic_faces (dof_handler, b_id, direction, offset); + + std::map 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 void identify_periodic_face_pairs @@ -2552,72 +2680,26 @@ next_cell: 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, - Point > 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::faces_per_cell; ++f) - if(cell->face(f)->boundary_indicator() == b_id1) - { - Point face_center (cell->face(f)->center()); - face_center(direction)=0.; - const std::pair - 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::faces_per_cell;++f) - if(cell->face(f)->boundary_indicator() == b_id2) - { - typename std::map, - Point >::iterator p - = face_locations.begin(); - - Point 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 - 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 > + 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 + periodic_tuple (it->cell[0], it->face_idx[0], + it->cell[1], it->face_idx[1]); + + periodicity_vector.push_back(periodic_tuple); + } } diff --git a/deal.II/source/grid/grid_tools.inst.in b/deal.II/source/grid/grid_tools.inst.in index fdaecb93c5..465931c953 100644 --- a/deal.II/source/grid/grid_tools.inst.in +++ b/deal.II/source/grid/grid_tools.inst.in @@ -234,6 +234,23 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II #if deal_II_dimension >= 2 + template + std::vector > + collect_periodic_faces + (const X &, + const types::boundary_id, + const types::boundary_id, + unsigned int, + const Tensor<1,deal_II_space_dimension> &); + + template + std::vector > + collect_periodic_faces + (const X &, + const types::boundary_id, + unsigned int, + const Tensor<1,deal_II_space_dimension> &); + template std::map > > collect_periodic_face_pairs @@ -243,7 +260,6 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II const types::boundary_id, int, const Tensor<1,deal_II_space_dimension> &); - template std::map > > collect_periodic_face_pairs @@ -270,7 +286,7 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II const types::boundary_id, int, const Tensor<1,deal_II_space_dimension> &); - + template void identify_periodic_face_pairs diff --git a/deal.II/source/hp/dof_handler.cc b/deal.II/source/hp/dof_handler.cc index d6afadf43c..fdb3a42ae5 100644 --- a/deal.II/source/hp/dof_handler.cc +++ b/deal.II/source/hp/dof_handler.cc @@ -1458,7 +1458,7 @@ namespace internal // 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; linen_raw_lines(); ++line) -- 2.39.5