Vector<double> &cell_vector);
void assemble_Neumann_boundary_terms(const FEFaceValues<dim> &face_fe,
- FullMatrix<double> &local_matrix,
- Vector<double> &local_vector);
+ FullMatrix<double> &local_matrix,
+ Vector<double> &local_vector);
void assemble_Dirichlet_boundary_terms(const FEFaceValues<dim> &face_fe,
- FullMatrix<double> &local_matrix,
- Vector<double> &local_vector,
- const double & h);
-
- void assemble_flux_terms(const FEFaceValuesBase<dim> &fe_face_values,
- const FEFaceValuesBase<dim> &fe_neighbor_face_values,
- FullMatrix<double> &vi_ui_matrix,
- FullMatrix<double> &vi_ue_matrix,
- FullMatrix<double> &ve_ui_matrix,
- FullMatrix<double> &ve_ue_matrix,
- const double & h);
+ FullMatrix<double> &local_matrix,
+ Vector<double> &local_vector,
+ const double & h);
+
+ void assemble_flux_terms(const FEFaceValuesBase<dim> &fe_face_values,
+ const FEFaceValuesBase<dim> &fe_neighbor_face_values,
+ FullMatrix<double> &vi_ui_matrix,
+ FullMatrix<double> &vi_ue_matrix,
+ FullMatrix<double> &ve_ui_matrix,
+ FullMatrix<double> &ve_ue_matrix,
+ const double & h);
void distribute_local_flux_to_global(
- const FullMatrix<double> & vi_ui_matrix,
- const FullMatrix<double> & vi_ue_matrix,
- const FullMatrix<double> & ve_ui_matrix,
- const FullMatrix<double> & ve_ue_matrix,
- const std::vector<types::global_dof_index> & local_dof_indices,
- const std::vector<types::global_dof_index> & local_neighbor_dof_indices);
+ const FullMatrix<double> & vi_ui_matrix,
+ const FullMatrix<double> & vi_ue_matrix,
+ const FullMatrix<double> & ve_ui_matrix,
+ const FullMatrix<double> & ve_ue_matrix,
+ const std::vector<types::global_dof_index> & local_dof_indices,
+ const std::vector<types::global_dof_index> & local_neighbor_dof_indices);
void solve();
// The constructor and destructor for this class is very much like the
// like those for step-40. The difference being that we'll be passing
// in an integer, <code>degree</code>, which tells us the maxiumum order
-// of the polynomial to use as well as <code>n_refine</code> which is the global
-// number of times we refine our mesh. The other main differences are
-// that we use a FESystem object for our choice of basis
+// of the polynomial to use as well as <code>n_refine</code> which is the
+// global number of times we refine our mesh. The other main differences
+// are that we use a FESystem object for our choice of basis
// functions. This is reminiscent of the mixed finite element method in
// step-20, however, in our case we use a FESystem
// of the form,
{
if((cell->center()[1]) > 0.9 )
{
- if((cell->center()[0] > 0.9) || (cell->center()[0] < 0.1) )
+ if((cell->center()[0] > 0.9) || (cell->center()[0] < 0.1))
cell->set_refine_flag();
}
}
// in order handle hanging nodes is no longer necessary. However,
// we will continue to use the constraint matrices inorder to efficiently
// distribute local computations to the global system, i.e. to the
- // <code>system_matrix</code> and <code>system_rhs</code>. Therefore, we just
- // instantiate the constraints matrix object, clear and close it.
+ // <code>system_matrix</code> and <code>system_rhs</code>. Therefore, we
+ // just instantiate the constraints matrix object, clear and close it.
constraints.clear();
constraints.close();
const unsigned int dofs_per_cell = fe.dofs_per_cell;
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
- std::vector<types::global_dof_index> local_neighbor_dof_indices(dofs_per_cell);
+ std::vector<types::global_dof_index>
+ local_neighbor_dof_indices(dofs_per_cell);
// We first remark that we have the FEValues objects for
// the values of our cell basis functions as was done in most
// FEFaceValues object, <code>fe_face_values</code>,
// for evaluating the basis functions
// on one side of an element face as well as another FEFaceValues object,
- // <code>fe_neighbor_face_values</code>, for evaluting the basis functions on the
- // opposite side of the face, i.e. on the neighoring element's face. In
- // addition, we also introduce a FESubfaceValues object,
+ // <code>fe_neighbor_face_values</code>, for evaluting the basis functions
+ // on the opposite side of the face, i.e. on the neighoring element's face.
+ // In addition, we also introduce a FESubfaceValues object,
// <code>fe_subface_values</code>, that
- // will be used for dealing with faces that have multiple refinement levels,
- // i.e. hanging nodes. When we have to evaulate the fluxes across
+ // will be used for dealing with faces that have multiple refinement
+ // levels, i.e. hanging nodes. When we have to evaulate the fluxes across
// a face that multiple refinement levels, we need to evaluate the
// fluxes across all its childrens' faces; we'll explain this more when
// the time comes.
FEFaceValues<dim> fe_face_values(fe,face_quadrature_formula,
face_update_flags);
- FEFaceValues<dim> fe_neighbor_face_values(fe, face_quadrature_formula,
- face_update_flags);
+ FEFaceValues<dim> fe_neighbor_face_values(fe,
+ face_quadrature_formula,
+ face_update_flags);
FESubfaceValues<dim> fe_subface_values(fe, face_quadrature_formula,
face_update_flags);
// and construct the local contribtuions from the numerical fluxes.
// The numerical fluxes will be due to 3 contributions: the
// interior faces, the faces on the Neumann boundary and the faces
- // on the Dirichlet boundary. We instantate a <code>face_iterator</code>
- // to loop
+ // on the Dirichlet boundary. We instantate a
+ // <code>face_iterator</code> to loop
// over all the faces of this cell and first see if the face is on
// the boundary. Notice how we do not reinitiaize the
// <code>fe_face_values</code>
- // object for the face until we know that we are actually on face that
- // lies on the boundary of the domain. The reason for doing this is for
- // computational efficiency; reinitializing the FEFaceValues for each face
- // is expensive and we do not want to do it unless we are actually
- // going use it to do computations. After this, we test if the face
+ // object for the face until we know that we are actually on face
+ // that lies on the boundary of the domain. The reason for doing this
+ // is for computational efficiency; reinitializing the FEFaceValues
+ // for each face is expensive and we do not want to do it unless we
+ // are actually going use it to do computations. After this, we test
+ // if the face
// is on the a Dirichlet or a Neumann segment of the boundary and
// call the appropriate subroutine to assemble the contributions for
// that boundary. Note that this assembles the flux contribution
face_no< GeometryInfo<dim>::faces_per_cell;
face_no++)
{
- typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+ typename DoFHandler<dim>::face_iterator face =
+ cell->face(face_no);
if(face->at_boundary() )
{
}
else
{
- // At this point we know that the face we are on is an interior
- // face. We can begin to assemble the interior flux matrices, but
- // first we want to make sure that the neighbor cell to this face
- // is a valid cell. Once we know that the neighbor is a valid cell
- // then we also want to get the meighbor cell that shares this cell's
- // face.
- Assert(cell->neighbor(face_no).state() == IteratorState::valid,
- ExcInternalError());
+ // At this point we know that the face we are on is an
+ // interior face. We can begin to assemble the interior
+ // flux matrices, but first we want to make sure that the
+ // neighbor cell to this face is a valid cell. Once we know
+ // that the neighbor is a valid cell then we also want to get
+ // the meighbor cell that shares this cell's face.
+ //
+ Assert(cell->neighbor(face_no).state() ==
+ IteratorState::valid,
+ ExcInternalError());
typename DoFHandler<dim>::cell_iterator neighbor =
cell->neighbor(face_no);
- // Now that we have the two cells whose face we want to compute
- // the numerical flux across, we need to know if the face has
- // been refined, i.e. if it has children faces. This occurs when
- // one of the cells has a different level of refinement than
- // the other cell. If this is the case, then this face has a
- // different level of refinement than the other faces of the cell,
- // i.e. on this face there is a hanging node. Hanging nodes are
- // not a problem in DG methods, the only time we have to watch
- // out for them is at this step and as you will see the changes
- // we have to our make are minor.
+ // Now that we have the two cells whose face we want to
+ // compute the numerical flux across, we need to know
+ // if the face has been refined, i.e. if it has children
+ // faces. This occurs when one of the cells has a
+ // different level of refinement than
+ // the other cell. If this is the case, then this face
+ // has a different level of refinement than the other faces
+ // of the cell, i.e. on this face there is a hanging node.
+ // Hanging nodes are not a problem in DG methods, the only
+ // time we have to watch out for them is at this step
+ // and as you will see the changes we have to our make
+ // are minor.
if(face->has_children())
{
- // We now need to find the face of our neighbor cell such that
- // neighbor(neigh_face_no) = cell(face_no).
+ // We now need to find the face of our neighbor cell
+ // such that neighbor(neigh_face_no) = cell(face_no).
const unsigned int neighbor_face_no =
cell->neighbor_of_neighbor(face_no);
- // Once we do this we then have to loop over all the subfaces
- // (children faces) of our cell's face and compute the interior
- // fluxes across the children faces and the neighbor's face.
+ // Once we do this we then have to loop over all the
+ // subfaces (children faces) of our cell's face and
+ // compute the interior fluxes across the children faces
+ // and the neighbor's face.
for(unsigned int subface_no=0;
subface_no < face->number_of_children();
++subface_no)
{
- // We then get the neighbor cell's subface that matches our
- // cell face's subface and the specific subface number.
- // We assert that the parent face cannot be more than one
- // level of refinement above the child's face. This is
- // because the deal.ii library does not allow neighboring
- // cells to have refinement levels that are more than
- // one level in difference.
+ // We then get the neighbor cell's subface that
+ // matches our cell face's subface and the
+ // specific subface number. We assert that the parent
+ // face cannot be more than one Level of
+ // refinement above the child's face. This is
+ // because the deal.ii library does not allow
+ // neighboring cells to have refinement levels
+ // that are more than one level in difference.
typename DoFHandler<dim>::cell_iterator neighbor_child =
- cell->neighbor_child_on_subface(face_no, subface_no);
+ cell->neighbor_child_on_subface(face_no,
+ subface_no);
- Assert(!neighbor_child->has_children(), ExcInternalError());
+ Assert(!neighbor_child->has_children(),
+ ExcInternalError());
- // Now that we are ready to build the local flux matrices for
+ // Now that we are ready to build the local flux
+ // matrices for
// this face we reset them e zero and
// reinitialize this <code>fe_values</code>
// to this cell's subface and
fe_neighbor_face_values.reinit(neighbor_child,
neighbor_face_no);
- // In addition, we get the minimum of diameters of the two cells
- // to include in the penalty term.
+ // In addition, we get the minimum of diameters of
+ // the two cells to include in the penalty term
double h = std::min(cell->diameter(),
neighbor_child->diameter());
- // We now finally assemble the interior fluxes for the
- // case of a face which has been refined using exactly
- // the same subroutine as we do when both cells have
- // the same refinement level.
+ // We now finally assemble the interior fluxes for
+ // the case of a face which has been refined using
+ // exactly the same subroutine as we do when both
+ // cells have the same refinement level.
assemble_flux_terms(fe_subface_values,
fe_neighbor_face_values,
vi_ui_matrix,
ve_ue_matrix,
h);
- // Now all that is left to be done before distribuing the
- // local flux matrices to the global system is get the neighbor
- // child faces dof indices.
+ // Now all that is left to be done before distribuing
+ // the local flux matrices to the global system
+ // is get the neighbor child faces dof indices.
neighbor_child->get_dof_indices(local_neighbor_dof_indices);
- // Once we have this cells dof indices and the neighboring
- // cell's dof indices we can use the ConstraintMatrix to
- // distribute the local flux matrices to the global system
- // matrix. This is done through the class function
+ // Once we have this cells dof indices and the
+ // neighboring cell's dof indices we can use the
+ // ConstraintMatrix to distribute the local flux
+ // matrices to the global system matrix.
+ // This is done through the class function
// <code>distribute_local_flux_to_global()</code>.
- distribute_local_flux_to_global(vi_ui_matrix,
- vi_ue_matrix,
- ve_ui_matrix,
- ve_ue_matrix,
- local_dof_indices,
- local_neighbor_dof_indices);
+ distribute_local_flux_to_global(
+ vi_ui_matrix,
+ vi_ue_matrix,
+ ve_ui_matrix,
+ ve_ue_matrix,
+ local_dof_indices,
+ local_neighbor_dof_indices);
}
}
else
{
- // At this point we know that this cell and the neighbor of
- // this cell are on the same refinement level and the work
- // to assemble the interior flux matrices is very much the
- // same as before. Infact it is much simpler since we do
- // not have to loop through the subfaces. However, we do
- // have to check that we do not compute the same contribution
- // twice. Since we are
- // looping over all the faces of all the cells in the mesh,
+ // At this point we know that this cell and the neighbor
+ // of this cell are on the same refinement level and
+ // the work to assemble the interior flux matrices
+ // is very much the same as before. Infact it is
+ // much simpler since we do not have to loop through the
+ // subfaces. However, we do have to check that we do
+ // not compute the same contribution twice. Since we are
+ // looping over all the faces of all the cells in the mesh,
// we pass over each face twice. If we do not take this
// into consideration when assembling the interior flux
// matrices we might compute the local interior flux matrix
// twice. To avoid doing this we only compute the interior
- // fluxes once for each face by restricting that the following
- // computation only occur on the on the cell face with
- // the lower index number.
+ // fluxes once for each face by restricting that the
+ // following computation only occur on the on
+ // the cell face with the lower index number.
if(neighbor->level() == cell->level() &&
neighbor->index() > cell->index() )
{
// Here we find the neighbor face such that
- // neighbor(neigh_face_no) = cell(face_no). In addition
- // we, reinitialize the FEFaceValues and neighbor cell's
- // FEFaceValues on their respective cells' faces,
- // as well as get the minimum diameter of this cell
- // and the neighbor cell and assign it to <code>h</code>.
+ // neighbor(neigh_face_no) = cell(face_no).
+ // In addition we, reinitialize the FEFaceValues
+ // and neighbor cell's FEFaceValues on their
+ // respective cells' faces, as well as get the
+ // minimum diameter of this cell
+ // and the neighbor cell and assign
+ // it to <code>h</code>.
const unsigned int neighbor_face_no =
cell->neighbor_of_neighbor(face_no);
ve_ue_matrix = 0;
fe_face_values.reinit(cell, face_no);
- fe_neighbor_face_values.reinit(neighbor, neighbor_face_no);
+ fe_neighbor_face_values.reinit(neighbor,
+ neighbor_face_no);
double h = std::min(cell->diameter(),
neighbor->diameter());
- // Just as before we assemble the interior fluxes using the
+ // Just as before we assemble the interior fluxes
+ // using the
// <code>assemble_flux_terms</code> subroutine,
// get the neighbor cell's
// face dof indices and use the constraint matrix to
// distribute the local flux matrices to the global
- // <code>system_matrix</code> using the class function
+ // <code>system_matrix</code> using the class
+ // function
// <code>distribute_local_flux_to_global()</code>
assemble_flux_terms(fe_face_values,
fe_neighbor_face_values,
neighbor->get_dof_indices(local_neighbor_dof_indices);
- distribute_local_flux_to_global(vi_ui_matrix,
- vi_ue_matrix,
- ve_ui_matrix,
- ve_ue_matrix,
- local_dof_indices,
- local_neighbor_dof_indices);
+ distribute_local_flux_to_global(
+ vi_ui_matrix,
+ vi_ue_matrix,
+ ve_ui_matrix,
+ ve_ue_matrix,
+ local_dof_indices,
+ local_neighbor_dof_indices);
}
{
for(unsigned int i=0; i<dofs_per_cell; i++)
{
- const Tensor<1, dim> psi_i_field = cell_fe[VectorField].value(i,q);
- const double div_psi_i_field = cell_fe[VectorField].divergence(i,q);
- const double psi_i_potential = cell_fe[Potential].value(i,q);
- const Tensor<1, dim> grad_psi_i_potential = cell_fe[Potential].gradient(i,q);
+ const Tensor<1, dim> psi_i_field = cell_fe[VectorField].value(i,q);
+ const double div_psi_i_field = cell_fe[VectorField].divergence(i,q);
+ const double psi_i_potential = cell_fe[Potential].value(i,q);
+ const Tensor<1, dim> grad_psi_i_potential = cell_fe[Potential].gradient(i,q);
for(unsigned int j=0; j<dofs_per_cell; j++)
{
{
for(unsigned int i=0; i<dofs_per_cell; i++)
{
- const Tensor<1, dim> psi_i_field = face_fe[VectorField].value(i,q);
- const double psi_i_potential = face_fe[Potential].value(i,q);
+ const Tensor<1, dim> psi_i_field = face_fe[VectorField].value(i,q);
+ const double psi_i_potential = face_fe[Potential].value(i,q);
for(unsigned int j=0; j<dofs_per_cell; j++)
{
- const Tensor<1, dim> psi_j_field = face_fe[VectorField].value(j,q);
- const double psi_j_potential = face_fe[Potential].value(j,q);
+ const Tensor<1, dim> psi_j_field = face_fe[VectorField].value(j,q);
+ const double psi_j_potential = face_fe[Potential].value(j,q);
// We compute contribution for the flux $\widehat{q}$ on
// the Dirichlet boundary which enters our system matrix as,
+
(penalty/h) *
psi_j_potential) *
- face_fe.JxW(q);
+ face_fe.JxW(q);
}
{
for(unsigned int i=0; i<dofs_per_cell; i++)
{
- const Tensor<1, dim> psi_i_field = face_fe[VectorField].value(i,q);
- const double psi_i_potential = face_fe[Potential].value(i,q);
+ const Tensor<1, dim> psi_i_field = face_fe[VectorField].value(i,q);
+ const double psi_i_potential = face_fe[Potential].value(i,q);
for(unsigned int j=0; j<dofs_per_cell; j++)
{
- const double psi_j_potential = face_fe[Potential].value(j,q);
+ const double psi_j_potential = face_fe[Potential].value(j,q);
// We compute contribution for the flux $\widehat{u}$ on the
// Neumann boundary which enters our system matrix as,
// the interior as,
//
// $\int_{\text{face}}
- // \left( \frac{1}{2} \, n^{-} \cdot ( \textbf{w}^{-} u^{-}
+ // \left( \frac{1}{2} \, n^{-}
+ // \cdot ( \textbf{w}^{-} u^{-}
// + w^{-} \textbf{q}^{-})
// + \boldsymbol \beta \cdot \textbf{w}^{-} u^{-}
// - w^{-} \boldsymbol \beta \cdot \textbf{q}^{-}
// to the computation,
//
// $\int_{\text{face}}
- // \left( \frac{1}{2} \, n^{-} \cdot ( \textbf{w}^{-} u^{+}
+ // \left( \frac{1}{2} \, n^{-} \cdot
+ // ( \textbf{w}^{-} u^{+}
// + w^{-} \textbf{q}^{+})
// - \boldsymbol \beta \cdot \textbf{w}^{-} u^{+}
// + w^{-} \boldsymbol \beta \cdot \textbf{q}^{+}
// to the computation,
//
// $ \int_{\text{face}}
- // \left( -\frac{1}{2}\, n^{-} \cdot (\textbf{w}^{+} u^{-}
+ // \left( -\frac{1}{2}\, n^{-} \cdot
+ // (\textbf{w}^{+} u^{-}
// + w^{+} \textbf{q}^{-} )
// - \boldsymbol \beta \cdot \textbf{w}^{+} u^{-}
// + w^{+} \boldsymbol \beta \cdot \textbf{q}^{-}
// cell to this face. This corresponds to the computation,
//
// $\int_{\text{face}}
- // \left( -\frac{1}{2}\, n^{-} \cdot ( \textbf{w}^{+} u^{+}
+ // \left( -\frac{1}{2}\, n^{-} \cdot
+ // ( \textbf{w}^{+} u^{+}
// + w^{+} \textbf{q}^{+} )
// + \boldsymbol \beta \cdot \textbf{w}^{+} u^{+}
// - w^{+} \boldsymbol \beta \cdot \textbf{q}^{+}
void
LDGPoissonProblem<dim>::
distribute_local_flux_to_global(
- const FullMatrix<double> & vi_ui_matrix,
- const FullMatrix<double> & vi_ue_matrix,
- const FullMatrix<double> & ve_ui_matrix,
- const FullMatrix<double> & ve_ue_matrix,
- const std::vector<types::global_dof_index> & local_dof_indices,
- const std::vector<types::global_dof_index> & local_neighbor_dof_indices)
+ const FullMatrix<double> & vi_ui_matrix,
+ const FullMatrix<double> & vi_ue_matrix,
+ const FullMatrix<double> & ve_ui_matrix,
+ const FullMatrix<double> & ve_ue_matrix,
+ const std::vector<types::global_dof_index> & local_dof_indices,
+ const std::vector<types::global_dof_index> & local_neighbor_dof_indices)
{
constraints.distribute_local_to_global(vi_ui_matrix,
local_dof_indices,
const std::string filename = ("solution." +
Utilities::int_to_string(
- triangulation.locally_owned_subdomain(),4));
+ triangulation.locally_owned_subdomain(),4));
std::ofstream output((filename + ".vtu").c_str());
data_out.write_vtu(output);
# Distributed Local Discontinuous Galerkin Methods
## Introduction
-This code is designed to numerically solve the <a href="https://en.wikipedia.org/wiki/Poisson's_equation">Poisson equation</a>,
+This code is designed to numerically solve the
+<a href="https://en.wikipedia.org/wiki/Poisson's_equation">Poisson equation</a>
@f{align}
-- \nabla \cdot \left(\ \nabla u \ \right)&= f(\textbf{x}) && \mbox{in} \ \Omega, \nonumber \\
--\nabla u \cdot \textbf{n} &= g_{N}(\textbf{x}) && \mbox{on} \ \partial \Omega_{N} \nonumber\\
+- \nabla \cdot \left(\ \nabla u \ \right)&= f(\textbf{x}) && \mbox{in} \
+\Omega,\nonumber \\
+-\nabla u \cdot \textbf{n} &= g_{N}(\textbf{x}) && \mbox{on} \ \partial
+\Omega_{N} \nonumber\\
u &= g_{D}(\textbf{x}) && \mbox{on} \ \partial \Omega_{D}. \nonumber
@f}
-in 2D and 3D using the local discontinuous Galerkin (LDG) method</a> from scratch. The tutorial codes, step-12 and step-39 use the MeshWorker interface to build discontinuous Galerkin (DG) methods. While this is very convenient, I could not use this framework for solving my problem and I needed write the LDG method from scratch. I thought it would help for others to have access to this example that goes through writing a discontinuous Galerkin method from scatch and also shows how to do it in a distributed setting using the <a href="https://www.trilinos.org">Trilinos</a> library. This example may also be of interest to users that wish to use the LDG method, as the method is distinctly different from the Interior Penalty Discontinuous Galerkin (IPDG) methods and was not covered in other tutorials on DG methods. The LDG method is very useful when one is working with a differential equation and desires both approximations to the scalar unknown function as well as its flux. The mixed method is another method where one can obtain both the scalar unknown function as we as its flux, however, the LDG method has less degrees of freedom when using the Raviart-Thomas elements in the mixed finite element. It also approximates the scalar unknown function and its flux using discontinous polynomial basis functions and are much more suitable when one wishes to use local refinement.
+in 2D and 3D using the local discontinuous Galerkin (LDG) method</a> from
+scratch. The tutorial codes step-12 and step-39 use the MeshWorker interface
+to build discontinuous Galerkin (DG) methods. While this is very convenient,
+I could not use this framework for solving my research problem and I
+needed write the LDG method from scratch. I thought it
+would be helpful for others to have access to
+this example that goes through writing a discontinuous Galerkin method from
+scatch and also shows how to do it in a distributed setting using the
+<a href="https://www.trilinos.org">Trilinos</a> library. This example may also
+be of interest to users that wish to use the LDG method, as the method is
+distinctly different from the Interior Penalty Discontinuous Galerkin (IPDG)
+methods and was not covered in other tutorials on DG methods. The LDG method
+is very useful when one is working with a differential equation and desires
+both approximations to the scalar unknown function as well as its flux. The
+mixed method is another method where one can obtain both the scalar unknown
+function as well as its flux, however, the LDG method has less degrees of
+freedom compared to the mixed method with the Raviart-Thomas element. It also
+approximates the scalar unknown function and its flux using discontinuous
+polynomial basis functions and are much more suitable when one wishes to use
+local refinement.
## Compiling and Running
-To generate a makefile for this code using CMake, type the following command into the terminal from the main directory:
+To generate a makefile for this code using CMake, type the following command
+into the terminal from the main directory:
cmake . -DDEAL_II_DIR=/path/to/deal.II
make release
-Either of these commands will create the executable, <code>main</code>, however the release mode will make a faster executable.
+Either of these commands will create the executable, <code>main</code>,
+however the release mode will make a faster executable.
-To run the code on <code>N</code> processors type the following command into ther terminal from the main directory,
+To run the code on <code>N</code> processors type the following command into
+the terminal from the main directory,
mpirun -np N ./main
-The output of the code will be in <code>.vtu</code> and <code>.pvtu</code> format and be written to disk in parallel. The results can be viewed using <a href="http://www.paraview.org/">ParaView</a>.
+The output of the code will be in <code>.vtu</code> and <code>.pvtu</code>
+format and be written to disk in parallel. The results can be viewed using
+<a href="http://www.paraview.org/">ParaView</a>.
## Local Discontinuous Galerkin Method
-In this section we discuss the LDG method and first introduce some notation. Let $\mathcal{T}_{h} = \mathcal{T}_{h}(\Omega) \, = \, \left\{ \, \Omega_{e} \, \right\}_{e=1}^{N}$ be the general triangulation of a domain $\Omega \; \subset \; \mathbb{R}^{d}, \; d \, = \, 1, 2, 3$, into $N$ non-overlapping elements $\Omega_{e}$ of diameter $h_{e}$. The maximum size of the diameters of all the elements is $h = \max( \, h_{e}\, )$. We define $\mathcal{E}_{h}$ to be the set of all element faces and $\mathcal{E}_{h}^{i} $ to be the set of all interior faces of elements which do not intersect the total boundary $(\partial \Omega)$. We define $\mathcal{E}_{D}$ and $\mathcal{E}_{N}$ to be the sets of all element faces and on the Dirichlet and Neumann boundaries respectively. Let $\partial \Omega_{e} \in \mathcal{E}_{h}^{i}$ be a interior boundary face element, we define the unit normal vector to be,
+In this section we discuss the LDG method and first introduce some notation.
+Let $\mathcal{T}_{h} = \mathcal{T}_{h}(\Omega) \, = \, \left\{ \, \Omega_{e}
+\, \right\}_{e=1}^{N}$ be the general triangulation of a domain $\Omega \;
+\subset \; \mathbb{R}^{d}, \; d \, = \, 1, 2, 3$, into $N$ non-overlapping
+elements $\Omega_{e}$ of diameter $h_{e}$. The maximum size of the diameters
+of all the elements is $h = \max( \, h_{e}\, )$. We define $\mathcal{E}_{h}$
+to be the set of all element faces and $\mathcal{E}_{h}^{i} $ to be the set of
+all interior faces of elements which do not intersect the total boundary
+$(\partial \Omega)$. We define $\mathcal{E}_{D}$ and $\mathcal{E}_{N}$ to be
+the sets of all element faces and on the Dirichlet and Neumann boundaries
+respectively. Let $\partial \Omega_{e} \in \mathcal{E}_{h}^{i}$ be a interior
+boundary face element, we define the unit normal vector to be,
@f{align}
-\textbf{n} \; = \; \text{unit normal vector to } \partial \Omega_{e} \text{ pointing from } \Omega_{e}^{-} \, \rightarrow \, \Omega_{e}^{+}.
+\textbf{n} \; = \; \text{unit normal vector to } \partial \Omega_{e}
+\text{ pointing from } \Omega_{e}^{-} \, \rightarrow \, \Omega_{e}^{+}.
@f}
We take the following definition on limits of functions on element faces,
@f{align}
-w^{-} (\textbf{x} ) \, \vert_{\partial \Omega_{e} } \; = \; \lim_{s \rightarrow 0^{-}} \, w(\textbf{x} + s \textbf{n}), && w^{+} (\textbf{x} ) \, \vert_{\partial \Omega_{e} } \; = \; \lim_{s \rightarrow 0^{+}} \, w(\textbf{x} + s \textbf{n}).
+w^{-} (\textbf{x} ) \, \vert_{\partial \Omega_{e} } \; = \;
+\lim_{s \rightarrow 0^{-}} \, w(\textbf{x} + s \textbf{n}), &&
+w^{+} (\textbf{x} ) \, \vert_{\partial \Omega_{e} } \; = \;
+\lim_{s \rightarrow 0^{+}} \, w(\textbf{x} + s \textbf{n}).
@f}
We define the average and jump of a function across an element face as,
@f{align}
\{\textbf{f} \} \; = \; \frac{1}{2}(\textbf{f}^- + \textbf{f}^+),
\qquad \mbox{and}\qquad
-\left[ \textbf{f} \right] \; = \;\textbf{f}^+ \cdot \textbf{n}^+ + \textbf{f}^- \cdot \textbf{n}^- ,
+\left[ \textbf{f} \right] \; = \;\textbf{f}^+ \cdot \textbf{n}^+ +
+\textbf{f}^- \cdot \textbf{n}^- ,
@f}
-where $f$ is a scalar function and $\textbf{f}$ is vector-valued function. We note that for a faces that are on the boundary of the domain we have,
+where $f$ is a scalar function and $\textbf{f}$ is vector-valued function.
+We note that for a faces that are on the boundary of the domain we have,
@f{align}
\left[ f \right] \; = \; f \, \textbf{n}
\qquad \mbox{and}\qquad
@f}
-We denote the volume integrals and surface integrals using the $L^{2}(\Omega)$ inner products by $( \, \cdot \, , \, \cdot \, )_{\Omega}$ and $\langle \, \cdot \, , \, \cdot \, \rangle_{\partial \Omega}$ respectively.
+We denote the volume integrals and surface integrals using the $L^{2}(\Omega)$
+inner products by $( \, \cdot \, , \, \cdot \, )_{\Omega}$ and $\langle \,
+\cdot \, , \, \cdot \, \rangle_{\partial \Omega}$ respectively.
-As with the mixed finite element method, the LDG discretization requires the Poisson equations be written as a first-order system. We do this by introducing an auxiliary variable which we call the current flux variable $\textbf{q}$:
+As with the mixed finite element method, the LDG discretization requires the
+Poisson equations be written as a first-order system. We do this by
+introducing an auxiliary variable which we call the current flux variable
+$\textbf{q}$:
@f{align}
\nabla \cdot \textbf{q}
\; &= \;
u &= g_{D}(\textbf{x}) && \mbox{on}\ \partial \Omega_{D}.
@f}
-In our numerical methods we will use approximations to scalar valued functions that reside in the finite-dimensional broken Sobolev spaces,
+In our numerical methods we will use approximations to scalar valued functions
+that reside in the finite-dimensional broken Sobolev spaces,
@f{align}
W_{h,k}
\, &= \,
-\left\{ w \in L^{2}(\Omega) \, : \; w \vert_{\Omega_{e}} \in \mathcal{Q}_{k,k}(\Omega_{e}), \quad \forall \, \Omega_{e} \in \mathcal{T}_{h} \right\},
+\left\{ w \in L^{2}(\Omega) \, : \; w \vert_{\Omega_{e}} \in
+\mathcal{Q}_{k,k}(\Omega_{e}), \quad \forall \, \Omega_{e} \in \mathcal{T}_{h}
+\right\},
@f}
-where $\mathcal{Q}_{k,k}(\Omega_{e})$ denotes the tensor product of discontinuous polynomials of order $k$ on the element $\Omega_{e}$. We use approximations of vector valued functions that are in,
+where $\mathcal{Q}_{k,k}(\Omega_{e})$ denotes the tensor product of
+discontinuous polynomials of order $k$ on the element $\Omega_{e}$. We use
+approximations of vector valued functions that are in,
@f{align}
\textbf{W}_{h,k}
\, &= \,
\left\{ \textbf{w} \in \left(L^{2}(\Omega)\right)^{d} \, :
- \; \textbf{w} \vert_{\Omega_{e}} \in \left( \mathcal{Q}_{k,k}(\Omega_{e}) \right)^{d}, \quad \forall \, \Omega_{e} \in \mathcal{T}_{h} \right\}
+ \; \textbf{w} \vert_{\Omega_{e}} \in \left( \mathcal{Q}_{k,k}(\Omega_{e})
+ \right)^{d}, \quad \forall \, \Omega_{e} \in \mathcal{T}_{h} \right\}
@f}
-We seek approximations for densities $u_{h} \in W_{h,k}$ and gradients $\textbf{q}_{h}\in \textbf{W}_{h,k}$. Multiplying (6) by $w \in W_{h,k}$ and (7) by $\textbf{w} \in \textbf{W}_{h,k}$ and integrating the divergence terms by parts over an element $\Omega_{e} \in \mathcal{T}_{h}$ we obtain,
+We seek approximations for densities $u_{h} \in W_{h,k}$ and gradients
+$\textbf{q}_{h}\in \textbf{W}_{h,k}$. Multiplying (6) by $w \in W_{h,k}$ and
+(7) by $\textbf{w} \in \textbf{W}_{h,k}$ and integrating the divergence terms
+by parts over an element $\Omega_{e} \in \mathcal{T}_{h}$ we obtain,
@f{align}
-
\left( \nabla w \, , \, \textbf{q}_{h} \right)_{\Omega_{e}}
-
\sum_{e} \left( \nabla w, \, \textbf{q}_{h} \right)_{\Omega_{e}}
+
-\langle \llbracket \, w \, \rrbracket \, , \, \widehat{\textbf{q}_{h} } \rangle_{\mathcal{E}_{h}^{i} }
+\langle \llbracket \, w \, \rrbracket \, , \, \widehat{\textbf{q}_{h} }
+\rangle_{\mathcal{E}_{h}^{i} }
+
-\langle \llbracket \, w \, \rrbracket \, , \, \widehat{\textbf{q}_{h} } \rangle_{\mathcal{E}_{D} \cup \mathcal{E}_{N}} \ &= \
+\langle \llbracket \, w \, \rrbracket \, , \, \widehat{\textbf{q}_{h} }
+\rangle_{\mathcal{E}_{D} \cup \mathcal{E}_{N}} \ &= \
\sum_{e} \left( w , \, f \right)_{\Omega_{e}} \\
\sum_{e} \left( \textbf{w} \, , \, \textbf{q}_{h} \right)_{\Omega_{e}}
-
\langle \, \llbracket \, \textbf{w} \, \rrbracket \, , \, \widehat{u_{h}}
\rangle_{\mathcal{E}_{h}^{i}}
+
-\langle \llbracket \, \textbf{w} \, \rrbracket \, , \, \widehat{u_{h}} \rangle_{\mathcal{E}_{D} \cup \mathcal{E}_{N}}
+\langle \llbracket \, \textbf{w} \, \rrbracket \, , \, \widehat{u_{h}}
+\rangle_{\mathcal{E}_{D} \cup \mathcal{E}_{N}}
\ &= \
0
@f}
for all $(w,\textbf{w}) \in W_{h,k} \times \textbf{W}_{h,k}$.
-The terms $\widehat{\textbf{q}_{h}}$ and $\widehat{u_{h}}$ are the numerical fluxes. The numerical fluxes are introduced to ensure consistency, stability, and enforce the boundary conditions weakly. The flux $\widehat{u_{h}}$ is,
+The terms $\widehat{\textbf{q}_{h}}$ and $\widehat{u_{h}}$ are the numerical
+fluxes. The numerical fluxes are introduced to ensure consistency, stability,
+and enforce the boundary conditions weakly. The flux $\widehat{u_{h}}$ is,
@f{align}
\widehat{u_{h}} \; = \; \left\{
\begin{array}{cl}
-\left\{ u_{h} \right\} \ + \ \boldsymbol \beta \cdot \left[ u_{h} \rrbracket & \ \text{in} \ \mathcal{E}_{h}^{i} \\
+\left\{ u_{h} \right\} \ + \ \boldsymbol \beta \cdot \left[ u_{h} \rrbracket &
+\ \text{in} \ \mathcal{E}_{h}^{i} \\
u_{h} & \ \text{in} \ \mathcal{E}_{N}\\
g_{D}(\textbf{x}) & \ \text{in} \ \mathcal{E}_{D} \\
\end{array}
@f{align}
\widehat{\textbf{q}_{h}} \; = \; \left\{
\begin{array}{cl}
-\left\{ \textbf{q}_{h} \right\} \ - \ \left[ \textbf{q}_{h} \right] \, \boldsymbol \beta \ + \ \sigma \, \left[ \, u_{h} \,
+\left\{ \textbf{q}_{h} \right\} \ - \ \left[ \textbf{q}_{h} \right] \,
+\boldsymbol \beta \ + \ \sigma \, \left[ \, u_{h} \,
\right] & \ \text{in} \ \mathcal{E}_{h}^{i} \\
g_{N}(\textbf{x}) \, \textbf{n} \, & \ \text{in} \ \mathcal{E}_{N}\\
-\textbf{q}_{h} \ + \ \sigma \, \left(u_{h} - g_{D}(\textbf{x}) \right) \, \textbf{n} & \ \text{in} \ \mathcal{E}_{D} \\
+\textbf{q}_{h} \ + \ \sigma \, \left(u_{h} - g_{D}(\textbf{x}) \right) \,
+\textbf{n} & \ \text{in} \ \mathcal{E}_{D} \\
\end{array}
\right.
@f}
-The term $\boldsymbol \beta$ is a constant unit vector which does not lie parallel to any element face in $ \mathcal{E}_{h}^{i}$. For $\boldsymbol \beta = 0$, $\widehat{\textbf{q}_{h}}$ and $\widehat{u_{h}}$ are called the central or Brezzi et. al. fluxes. For $\boldsymbol \beta \neq 0$, $\widehat{\textbf{q}_{h}}$ and $\widehat{u_{h}}$ are called the LDG/alternating fluxes. The term $\sigma$ is the penalty parameter that is defined as,
+The term $\boldsymbol \beta$ is a constant unit vector which does not lie
+parallel to any element face in $ \mathcal{E}_{h}^{i}$. For
+$\boldsymbol \beta = 0$, $\widehat{\textbf{q}_{h}}$ and $\widehat{u_{h}}$
+are called the central or Brezzi et. al. fluxes. For
+$\boldsymbol \beta \neq 0$, $\widehat{\textbf{q}_{h}}$ and $\widehat{u_{h}}$
+are called the LDG/alternating fluxes. The term $\sigma$ is the penalty
+parameter that is defined as,
@f{align}
\sigma \; = \; \left\{
\begin{array}{cc}
-\tilde{\sigma} \, \min \left( h^{-1}_{e_{1}}, h^{-1}_{e_{2}} \right) & \textbf{x} \in \langle \Omega_{e_{1}}, \Omega_{e_{2}} \rangle \\
-\tilde{\sigma} \, h^{-1}_{e} & \textbf{x} \in \partial \Omega_{e} \cap \in \mathcal{E}_{D}
+\tilde{\sigma} \, \min \left( h^{-1}_{e_{1}}, h^{-1}_{e_{2}} \right) &
+\textbf{x} \in \langle \Omega_{e_{1}}, \Omega_{e_{2}} \rangle \\
+\tilde{\sigma} \, h^{-1}_{e} & \textbf{x} \in \partial \Omega_{e} \cap
+\in \mathcal{E}_{D}
\end{array}
\right.
\label{eq:Penalty}
with $\tilde{\sigma}$ being a positive constant.
-We can now substitute (16) and (17) into (14) and (15) to obtain the solution pair $(u_{h}, \textbf{q}_{h})$ to the LDG approximation to the Poisson equation given by:
+We can now substitute (16) and (17) into (14) and (15) to obtain the solution
+pair $(u_{h}, \textbf{q}_{h})$ to the LDG approximation to the Poisson
+equation given by:
Find $u_{h} \in W_{h,k}$ and $\textbf{q}_{h} \in \textbf{W}_{h,k}$ such that,
@f{align}
-a(\textbf{w}, \textbf{q}_{h}) \ + \ b^{T}(\textbf{w}, u_{h}) \ &= \ G(\textbf{w}) \nonumber \\
+a(\textbf{w}, \textbf{q}_{h}) \ + \ b^{T}(\textbf{w}, u_{h}) \ &= \
+G(\textbf{w}) \nonumber \\
b(w, \textbf{q}_{h}) \ + \ c(w, u_{h}) \ &= \ F(w)
\label{eq:LDG_bilinear}
@f}
-for all $(w, \textbf{w}) \in W_{h,k} \times \textbf{W}_{h,k}$. This leads to the linear system,
+for all $(w, \textbf{w}) \in W_{h,k} \times \textbf{W}_{h,k}$. This leads to
+the linear system,
@f{align}
\left[
\right]
@f}
-Where $\textbf{U}$ and $\textbf{Q}$ are the degrees of freedom vectors for $u_{h}$ and $\textbf{q}_{h}$ respectively. The terms $\textbf{G}$ and $\textbf{F}$ are the corresponding vectors to $G(\textbf{w})$ and $F(w)$ respectively. The matrix in for the LDG system is non-singular for any $\sigma > 0$.
+Where $\textbf{U}$ and $\textbf{Q}$ are the degrees of freedom vectors for
+$u_{h}$ and $\textbf{q}_{h}$ respectively. The terms $\textbf{G}$ and
+$\textbf{F}$ are the corresponding vectors to $G(\textbf{w})$ and $F(w)$
+respectively. The matrix in for the LDG system is non-singular for any
+$\sigma > 0$.
The bilinear forms in (19) and right hand functions are defined as,
\sum_{e} \left(\nabla w, \textbf{q}_{h} \right)_{\Omega_{e}}
+
\langle \left[ w \right],
-\left\{\textbf{q}_{h} \right\} - \left[ \textbf{q}_{h} \right] \boldsymbol \beta \rangle_{\mathcal{E}_{h}^{i}}
+\left\{\textbf{q}_{h} \right\} - \left[ \textbf{q}_{h} \right] \boldsymbol
+\beta \rangle_{\mathcal{E}_{h}^{i}}
+
\langle w, \textbf{n} \cdot \textbf{q}_{h} \rangle_{\mathcal{E}_{D}}\\
a(\textbf{w},\textbf{q}_{h}) \, &= \,
-
\sum_{e} \left(\nabla \cdot \textbf{w}, u_{h} \right)_{\Omega_{e}}
+
-\langle \left[ \textbf{w} \right], \left\{u_{h} \right\} + \boldsymbol \beta \cdot \left[ u_{h} \right] \rangle_{\mathcal{E}_h^{i} }
+\langle \left[ \textbf{w} \right], \left\{u_{h} \right\} + \boldsymbol
+\beta \cdot \left[ u_{h} \right] \rangle_{\mathcal{E}_h^{i} }
+
\langle w, u_{h} \rangle_{\mathcal{E}_{N} } \\
c(w,u_{h}) \, &= \,
-\langle \left[ w \right], \sigma \left[ u_{h} \right] \rangle_{\mathcal{E}_{h}^{i}}
+\langle \left[ w \right], \sigma \left[ u_{h} \right]
+\rangle_{\mathcal{E}_{h}^{i}}
+
\langle w, \sigma u_{h} \rangle_{\mathcal{E}_{D}} \\
-G(\textbf{w}) \ & = \ - \langle \textbf{w}, g_{D} \rangle_{\mathcal{E}_{D}}\\
-F(w) \ & = \ \sum_{e} (w,f)_{\Omega_{e}} - \langle w, g_{N} \rangle_{\mathcal{E}_{N}} + \langle w, \sigma g_{D} \rangle_{\mathcal{E}_{D}}
+G(\textbf{w}) \ & = \ - \langle \textbf{w}, g_{D}
+\rangle_{\mathcal{E}_{D}}\\
+F(w) \ & = \ \sum_{e} (w,f)_{\Omega_{e}} - \langle w, g_{N}
+\rangle_{\mathcal{E}_{N}} + \langle w, \sigma g_{D} \rangle_{\mathcal{E}_{D}}
@f}
-As discussed in step-20, we won't be assembling the bilinear terms explicitly, instead we will assemble all the solid integrals and fluxes at once. We note that in order to actually build the flux terms in our local flux matrices we will substitute in the definitions in the bilinear terms above.
+As discussed in step-20, we won't be assembling the bilinear terms explicitly,
+instead we will assemble all the solid integrals and fluxes at once. We note
+that in order to actually build the flux terms in our local flux matrices we
+will substitute in the definitions in the bilinear terms above.
# The Commented Code
\ No newline at end of file
projects / code-gallery.git / commitdiff