+<br>
+
+<i>
+This program was contributed by Martin Kronbichler and Scott Miller.
+</i>
+
<a name="Intro"></a>
<h1>Introduction</h1>
and similarly for edges and faces. As another example,
for the FE_DGP_Monomial basis, each
scalar solution component is represented by polynomials of degree $p$
-which yields $(1/dim!)*\prod_{i=1}^{dim}(k+i)$ degrees of freedom per
+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
to all of the degrees of freedom in the adjacent elements. The resulting
discrete equations yield very large linear systems very quickly, especially
To alleviate the computational cost of solving such large linear systems,
the hybridizable discontinuous Galerkin (HDG) methodology was introduced
by Cockburn and co-workers
-(N.C. Nguyen and J. Peraire:
+(see the references in the recent HDG overview article by
+ N.C. Nguyen and J. Peraire:
<i>Hybridizable discontinuous Galerkin methods for partial differential
equations in continuum mechanics</i>, Journal of Computational Physics,
2012, 231:18, 5955-5988.
We multiply these equations by the weight functions $\mathbf{v}, w$
and integrate by parts over every element $K$ to obtain:
@f{eqnarray*}
- (\kappa^{-1} \mathbf{q}, \mathbf{v})_K - (u, \nabla\cdot\mathbf{v})_K
- + \left<\hat{u}, \mathbf{v}\cdot\mathbf{n}\right>_{\partial K} &=& 0, \\
- - (\mathbf{c} u + \mathbf{q}, \nabla w)_K
- + \left<(\hat{\mathbf{c} u}+\hat{\mathbf{q}})\cdot\mathbf{n}, w\right>_{\partial K}
- &=& (f,w)_K.
+ (\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}
+ &=& (w,f)_K.
@f}
The terms decorated with a hat denote the numerical traces (also commonly referred
@f}
Eliminating $\hat{u}$ from the weak form in favor of $\lambda$, summing the elemental
-contributions across all elements in the triangulation, and enforcing the normal
-component of the numerical flux, we arrive at the final form of the problem:
-Find $(\mathbf{q}, u, \lambda)$ such that
-TODO: update this equation!
+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
+\mathcal{V}_h^p \times \mathcal{W}_h^p \times \mathcal{M}_h^p$ such that
@f{eqnarray*}
- (\kappa^{-1} \mathbf{q}, \mathbf{v})_{\mathcal{T}}
- - (u, \nabla\cdot\mathbf{v})_{\mathcal{T}}
- + \left<\lambda, \mathbf{v}\cdot\mathbf{n}\right>_{\partial\mathcal{T}}
+ (\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<g_D, \mathbf{v}\cdot\mathbf{n}\right>_{\partial\Omega_D},
+ - \left<\mathbf{v}\cdot\mathbf{n}, g_D\right>_{\partial\Omega_D},
+ \quad \forall \mathbf{v} \in \mathcal{V}_h^p,
\\
- - (\mathbf{c} u + \mathbf{q}, \nabla w)_{\mathcal{T}}
- + \left<(\hat{\mathbf{c}u} + \hat{\mathbf{q}})\cdot\mathbf{n}, w\right>_{\partial \mathcal{T}}
+ (w, \mathbf{c}\nabla u_h + \nabla \cdot \mathbf{q}_h)_{\mathcal{T}}
+ + \left<(w, \tau (u_h - \lambda_h)\right>_{\partial \mathcal{T}}
&=&
- (f, w)_{\mathcal{T}},
+ (w, f)_{\mathcal{T}},
+ \quad \forall w \in \mathcal{W}_h^p,
\\
- < \left[ \negthinspace \left[ (\hat{\mathbf{c}u} + \hat{\mathbf{q}})\cdot\mathbf{n}
- \right] \negthinspace \right], \mu>_{\partial \mathcal{T}}
+ < \mu, \mathbf{c} \lambda_h\cdot \mathbf{n}
+ + \mathbf{q}_h\cdot \mathbf{n}
+ + \tau (u_h - \lambda_h)>_{\partial \mathcal{T}}
&=&
- <g_N, \mu>_{\partial\Omega_N}
+ <\mu, g_N>_{\partial\Omega_N},
+ \quad \forall \mu \in \mathcal{M}_h^p.
@f}
-for all $(\mathbf{v}, w, \mu)$.
+
+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
+variable which has support on the codimension-1 surfaces (faces) of the mesh.
+In the weak form given above, we can note the following coupling patterns:
+<ol>
+ <li> The matrix $A$ consists of local-local coupling terms. These arise when the
+ local weighting functions $(\mathbf{v}, w)$ multiply the local solution terms
+ $(\mathbf{q}_h, u_h)$.
+ <li> The matrix $B$ represents the local-face coupling. These are the terms
+ with weighting functions $(\mathbf{v}, w)$ multiplying the skeleton variable
+ $\lambda_h$.
+ <li> The matrix $C$ represents the face-local coupling, which involves the
+ weighting function $\mu$ multiplying the local solutions $(\mathbf{q}_h, u_h)$.
+ <li> The matrix $D$ is the face-face coupling;
+ terms involve both $\mu$ and $\lambda_h$.
+</ol>
+
+<h4> Post-processing and super-convergence </h4>
<h3> Problem specific data </h3>
// the simulation.
#include <deal.II/numerics/data_out_faces.h>
+namespace Step51
+{
+
using namespace dealii;
// @sect3{Equation data}
* this->width);
}
-// @sect3{The Step51 HDG solver class}
+// @sect3{The HDG HDG solver class}
// The HDG solution procedure follows closely that of step-7. The major
// difference is the use of 3 different sets of <code>DoFHandler</code> and FE
// solutions from the skeleton values) and once for the postprocessing where
// we extract a solution that converges at higher order.
template <int dim>
-class Step51
+class HDG
{
public:
enum RefinementMode
global_refinement, adaptive_refinement
};
- Step51 (const unsigned int degree,
+ HDG (const unsigned int degree,
const RefinementMode refinement_mode);
void run ();
ConvergenceTable convergence_table;
};
-// @sect3{The Step51 class implementation}
+// @sect3{The HDG class implementation}
-// @sect4{Step51::Step51}
+// @sect4{HDG::HDG}
// The constructor is similar to those in other examples, with the
// exception of handling multiple <code>DoFHandler</code> and
// <code>FiniteElement</code> objects.
template <int dim>
-Step51<dim>::Step51 (const unsigned int degree,
+HDG<dim>::HDG (const unsigned int degree,
const RefinementMode refinement_mode) :
fe_local (FE_DGQ<dim>(degree), dim,
FE_DGQ<dim>(degree), 1),
{}
-// @sect4{Step51::PerTaskData}
+// @sect4{HDG::PerTaskData}
// First 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.
template <int dim>
-struct Step51<dim>::PerTaskData
+struct HDG<dim>::PerTaskData
{
FullMatrix<double> cell_matrix;
Vector<double> cell_vector;
}
};
-// @sect4{Step51::ScratchData}
+// @sect4{HDG::ScratchData}
// @p ScratchData contains persistent data for each thread within
// <code>WorkStream</code>. The <code>FEValues</code>,
// matrix, and vector objects should be familiar by now.
// a large number of zero terms on each cell, which would significantly
// slow the program.
template <int dim>
-struct Step51<dim>::ScratchData
+struct HDG<dim>::ScratchData
{
FEValues<dim> fe_values_local;
FEFaceValues<dim> fe_face_values_local;
};
-// @sect4{Step51::PostProcessScratchData}
+// @sect4{HDG::PostProcessScratchData}
// @p PostProcessScratchData contains the data used by <code>WorkStream</code>
// when post-processing the local solution $u^*$. It is similar, but much
// simpler, than @p ScratchData.
template <int dim>
-struct Step51<dim>::PostProcessScratchData
+struct HDG<dim>::PostProcessScratchData
{
FEValues<dim> fe_values_local;
FEValues<dim> fe_values;
};
-// @sect4{Step51::copy_local_to_global}
+// @sect4{HDG::copy_local_to_global}
// If we are in the first step of the solution, i.e. @p trace_reconstruct=false,
// then we assemble the global system.
template <int dim>
-void Step51<dim>::copy_local_to_global(const PerTaskData &data)
+void HDG<dim>::copy_local_to_global(const PerTaskData &data)
{
if (data.trace_reconstruct == false)
constraints.distribute_local_to_global (data.cell_matrix,
system_matrix, system_rhs);
}
-// @sect4{Step51::setup_system}
+// @sect4{HDG::setup_system}
// The system for an HDG solution is setup in an analogous manner to most
// of the other tutorial programs. We are careful to distribute dofs with
// all of our <code>DoFHandler</code> objects. The @p solution and @p system_matrix
// objects go with the global skeleton solution.
template <int dim>
void
-Step51<dim>::setup_system ()
+HDG<dim>::setup_system ()
{
dof_handler_local.distribute_dofs(fe_local);
dof_handler.distribute_dofs(fe);
}
-// @sect4{Step51::assemble_system}
+// @sect4{HDG::assemble_system}
// The @p assemble_system function is similar to <code>Step-32</code>, where a few
// objects are setup, and then <code>WorkStream</code> is used to do the work in a
// multi-threaded manner. The @p trace_reconstruct input parameter is used
// global skeleton solution (false).
template <int dim>
void
-Step51<dim>::assemble_system (const bool trace_reconstruct)
+HDG<dim>::assemble_system (const bool trace_reconstruct)
{
const QGauss<dim> quadrature_formula(fe.degree+1);
const QGauss<dim-1> face_quadrature_formula(fe.degree+1);
WorkStream::run(dof_handler.begin_active(),
dof_handler.end(),
*this,
- &Step51<dim>::assemble_system_one_cell,
- &Step51<dim>::copy_local_to_global,
+ &HDG<dim>::assemble_system_one_cell,
+ &HDG<dim>::copy_local_to_global,
scratch,
task_data);
}
-// @sect4{Step51::assemble_system_one_cell}
-// The real work of the Step51 program is done by @p assemble_system_one_cell.
+// @sect4{HDG::assemble_system_one_cell}
+// The real work of the HDG program is done by @p assemble_system_one_cell.
// Assembling the local matrices $A, B, C$ is done here, along with the
// local contributions of the global matrix $D$.
template <int dim>
void
-Step51<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+HDG<dim>::assemble_system_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData &scratch,
PerTaskData &task_data)
{
tau_stab) * scratch.u_phi[i])
* scratch.tr_phi[j]
) * JxW;
+
+// Note the sign of the face-local matrix. We negate the sign during
+// assembly here so that we can use the FullMatrix::mmult with addition
+// when computing the Schur complement.
scratch.fl_matrix(jj,ii) -= (
(scratch.q_phi[i] * normal
+
}
-// @sect4{Step51::solve}
+// @sect4{HDG::solve}
// The skeleton solution is solved for by using a BiCGStab solver with
// identity preconditioner.
template <int dim>
-void Step51<dim>::solve ()
+void HDG<dim>::solve ()
{
SolverControl solver_control (system_matrix.m()*10,
1e-11*system_rhs.l2_norm());
template <int dim>
void
-Step51<dim>::postprocess()
+HDG<dim>::postprocess()
{
// construct post-processed solution with (hopefully) higher order of
// accuracy
WorkStream::run(dof_handler_u_post.begin_active(),
dof_handler_u_post.end(),
- std_cxx1x::bind (&Step51<dim>::postprocess_one_cell,
+ std_cxx1x::bind (&HDG<dim>::postprocess_one_cell,
std_cxx1x::ref(*this),
std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3),
std_cxx1x::function<void(const unsigned int&)>(),
template <int dim>
void
-Step51<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+HDG<dim>::postprocess_one_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
PostProcessScratchData &scratch,
unsigned int &)
{
cell->distribute_local_to_global(scratch.cell_sol, solution_u_post);
}
-// @sect4{Step51::output_results}
+// @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
// objects. The graphical output for the skeleton variable is done through
// use of the <code>DataOutFaces</code> class.
template <int dim>
-void Step51<dim>::output_results (const unsigned int cycle)
+void HDG<dim>::output_results (const unsigned int cycle)
{
std::string filename;
switch (refinement_mode)
data_out_face.write_vtk (face_output);
}
-// @sect4{Step51::refine_grid}
-// We implement two different refinement cases for Step51, just as in
+// @sect4{HDG::refine_grid}
+// We implement two different refinement cases for HDG, just as in
// <code>Step-7</code>: adaptive_refinement and global_refinement.
// The global_refinement option recreates the entire triangulation
// every time.
// The adaptive_refinement mode uses the <code>KellyErrorEstimator</code>
// give a decent approximation of the non-regularity of the local solutions.
template <int dim>
-void Step51<dim>::refine_grid (const unsigned int cycle)
+void HDG<dim>::refine_grid (const unsigned int cycle)
{
if (cycle == 0)
{
cell->face(face)->set_boundary_indicator (1);
}
-// @sect4{Step51::run}
+// @sect4{HDG::run}
// The functionality here is basically the same as <code>Step-7</code>.
// We loop over 10 cycles, refining the grid on each one. At the end,
// convergence tables are created.
template <int dim>
-void Step51<dim>::run ()
+void HDG<dim>::run ()
{
for (unsigned int cycle=0; cycle<10; ++cycle)
{
convergence_table.write_text(std::cout);
}
+}//Step51
int main (int argc, char **argv)
{
<< "=============================================" << std::endl
<< std::endl;
- Step51<dim> hdg_problem (1, Step51<dim>::adaptive_refinement);
+ Step51::HDG<dim> hdg_problem (1, Step51::HDG<dim>::adaptive_refinement);
hdg_problem.run ();
std::cout << std::endl;
<< "===========================================" << std::endl
<< std::endl;
- Step51<dim> hdg_problem (1, Step51<dim>::global_refinement);
+ Step51::HDG<dim> hdg_problem (1, Step51::HDG<dim>::global_refinement);
hdg_problem.run ();
std::cout << std::endl;
<< "===========================================" << std::endl
<< std::endl;
- Step51<dim> hdg_problem (3, Step51<dim>::global_refinement);
+ Step51::HDG<dim> hdg_problem (3, Step51::HDG<dim>::global_refinement);
hdg_problem.run ();
std::cout << std::endl;
projects / dealii-svn.git / commitdiff