--- /dev/null
+// @sect3{LDGPoisson.cc}
+// The code begins as per usual with a long list of the the included
+// files from the deal.ii library.
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+// Here's where the classes for the DG methods begin.
+// We can use either the Lagrange polynomials,
+#include <deal.II/fe/fe_dgq.h>
+// or the Legendre polynomials
+#include <deal.II/fe/fe_dgp.h>
+// as basis functions. I'll be using the Lagrange polynomials.
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+
+// Now we have to load in the deal.II files that will allow us to use
+// a distributed computing framework.
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/distributed/tria.h>
+
+
+// Additionally we load the files that will allow us to interact with
+// the Trilinos library.
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/trilinos_solver.h>
+
+
+// The functions class contains all the defintions of the functions we
+// will use, i.e. the right hand side function, the boundary conditions
+// and the test functions.
+#include "Functions.cc"
+
+using namespace dealii;
+
+
+// Here is the main class for the Local Discontinuous Galerkin method
+// applied to Poisson's equation, we won't explain much of the
+// the class and method declarations, but dive deeper into describing the
+// functions when the are defined. The only thing I will menion
+// about the class declaration is that here is where I labeled
+// the different types of boundaries using enums.
+template <int dim>
+class LDGPoissonProblem
+{
+
+public:
+ LDGPoissonProblem(const unsigned int degree,
+ const unsigned int n_refine);
+
+ ~LDGPoissonProblem();
+
+ void run();
+
+
+private:
+ void make_grid();
+
+ void make_dofs();
+
+ void assemble_system();
+
+ void assemble_cell_terms(const FEValues<dim> &cell_fe,
+ FullMatrix<double> &cell_matrix,
+ Vector<double> &cell_vector);
+
+ void assemble_Neumann_boundary_terms(const FEFaceValues<dim> &face_fe,
+ 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);
+
+ 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);
+
+ void solve();
+
+ void output_results() const;
+
+ const unsigned int degree;
+ const unsigned int n_refine;
+ double penalty;
+ double h_max;
+ double h_min;
+
+ enum
+ {
+ Dirichlet,
+ Neumann
+ };
+
+ parallel::distributed::Triangulation<dim> triangulation;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
+
+ ConstraintMatrix constraints;
+
+ SparsityPattern sparsity_pattern;
+
+ TrilinosWrappers::SparseMatrix system_matrix;
+ TrilinosWrappers::MPI::Vector locally_relevant_solution;
+ TrilinosWrappers::MPI::Vector system_rhs;
+
+ ConditionalOStream pcout;
+ TimerOutput computing_timer;
+
+ SolverControl solver_control;
+ TrilinosWrappers::SolverDirect solver;
+
+ const RightHandSide<dim> rhs_function;
+ const DirichletBoundaryValues<dim> Dirichlet_bc_function;
+ const TrueSolution<dim> true_solution;
+};
+
+
+// @sect4{Class constructor and destructor}
+// 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
+// functions. This is reminiscent of the mixed finite element method in
+// step-20, however, in our case we use a FESystem
+// of the form,
+//
+// <code>
+// fe( FESystem<dim>(FE_DGQ<dim>(degree), dim), 1,
+// FE_DGQ<dim>(degree), 1)
+// </code>
+//
+// which tells us that the basis functions contain discontinous polynomials
+// of order <code>degree</code> in each of the <code>dim</code> dimensions
+// for the vector field. For the scalar unknown we
+// use a discontinuous polynomial of the order <code>degree</code>.
+// The LDG method solves for both the primary variable as well as
+// its gradient, just like the mixed finite element method. However,
+// unlike the mixed method, the LDG method ues discontinuous
+// polynomials to approximate both variables.
+// The other difference bewteen our constructor and that of step-40 is that
+// we all instantiate our linear solver in the constructor definition.
+template <int dim>
+LDGPoissonProblem<dim>::
+LDGPoissonProblem(const unsigned int degree,
+ const unsigned int n_refine)
+ :
+ degree(degree),
+ n_refine(n_refine),
+ triangulation(MPI_COMM_WORLD,
+ typename Triangulation<dim>::MeshSmoothing
+ (Triangulation<dim>::smoothing_on_refinement |
+ Triangulation<dim>::smoothing_on_coarsening)),
+ fe( FESystem<dim>(FE_DGQ<dim>(degree), dim), 1,
+ FE_DGQ<dim>(degree), 1),
+ dof_handler(triangulation),
+ pcout(std::cout,
+ Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0),
+ computing_timer(MPI_COMM_WORLD,
+ pcout,
+ TimerOutput::summary,
+ TimerOutput::wall_times),
+ solver_control(1),
+ solver(solver_control),
+ rhs_function(),
+ Dirichlet_bc_function()
+{
+}
+
+template <int dim>
+LDGPoissonProblem<dim>::
+~LDGPoissonProblem()
+{
+ dof_handler.clear();
+}
+
+// @sect4{Make_grid}
+// This function shows how to make a grid using local
+// refinement and also shows how to label the boundaries
+// using the defined enum.
+template <int dim>
+void
+LDGPoissonProblem<dim>::
+make_grid()
+{
+ GridGenerator::hyper_cube(triangulation, 0, 1);
+ triangulation.refine_global(n_refine);
+
+ unsigned int local_refine = 2;
+ for(unsigned int i =0; i <local_refine; i++)
+ {
+ typename Triangulation<dim>::active_cell_iterator
+ cell = triangulation.begin_active(),
+ endc = triangulation.end();
+
+ // We loop over all the cells in the mesh
+ // and mark the appropriate cells for refinement.
+ // In this example we only choose cells which are
+ // near $x=0$ and $x=1$ in the
+ // the domain. This was just to show that
+ // the LDG method is working with local
+ // refinement and discussions on building
+ // more realistic refinement stategies are
+ // discussed elsewhere in the deal.ii
+ // documentation.
+ for(; cell != endc; cell++)
+ {
+ if((cell->center()[1]) > 0.9 )
+ {
+ if((cell->center()[0] > 0.9) || (cell->center()[0] < 0.1) )
+ cell->set_refine_flag();
+ }
+ }
+ // Now that we have marked all the cells
+ // that we want to refine locally we can go ahead and
+ // refine them.
+ triangulation.execute_coarsening_and_refinement();
+ }
+
+ // To label the boundary faces of the mesh with their
+ // type, i.e. Dirichlet or Neumann,
+ // we loop over all the cells in the mesh and then over
+ // all the faces of each cell. We then have to figure out
+ // which faces are on the bounadry and set all faces
+ // on the boundary to have
+ // <code>boundary_id</code> to be <code>Dirichlet</code>.
+ // We remark that one could easily set more complicated
+ // conditions where there are both Dirichlet or
+ // Neumann boundaries.
+ typename Triangulation<dim>::cell_iterator
+ cell = triangulation.begin(),
+ endc = triangulation.end();
+ for(; cell != endc; cell++)
+ {
+ for(unsigned int face_no=0;
+ face_no < GeometryInfo<dim>::faces_per_cell;
+ face_no++)
+ {
+ if(cell->face(face_no)->at_boundary() )
+ cell->face(face_no)->set_boundary_id(Dirichlet);
+ }
+ }
+}
+
+// @sect3{make_dofs}
+// This function is responsible for distributing the degrees of
+// freedom (dofs) to the processors and allocating memory for
+// the global system matrix, <code>system_matrix</code>,
+// and global right hand side vector, <code>system_rhs</code> .
+// The dofs are the unknown coefficients for the polynomial
+// approximation of our solution to Poisson's equation in the scalar
+// variable and its gradient.
+template <int dim>
+void
+LDGPoissonProblem<dim>::
+make_dofs()
+{
+ TimerOutput::Scope t(computing_timer, "setup");
+
+ // The first step to setting up our linear system is to
+ // distribute the degrees of freedom (dofs) across the processors,
+ // this is done with the <code>distribute_dofs()</code>
+ // method of the DoFHandler. We remark
+ // the same exact function call that occurs when using deal.ii
+ // on a single machine, the DoFHandler automatically knows
+ // that we are distributed setting because it was instantiated
+ // with a distributed triangulation!
+ dof_handler.distribute_dofs(fe);
+
+ // We now renumber the dofs so that the vector of unkonwn dofs
+ // that we are solving for, <code>locally_relevant_solution</code>
+ // corresponds to a vector of the form,
+ //
+ // $ \left[\begin{matrix} \textbf{Q} \\ \textbf{U} \end{matrix}\right] $
+ DoFRenumbering::component_wise(dof_handler);
+
+ // Now we get the locally owned dofs, that is the dofs that our local
+ // to this processor. These dofs corresponding entries in the
+ // matrix and vectors that we will write to.
+ IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs();
+
+ // In additon to the locally owned dofs, we also need the the locally
+ // relevant dofs. These are the dofs that have read access to and we
+ // need in order to do computations on our processor, but, that
+ // we do not have the ability to write to.
+ IndexSet locally_relevant_dofs;
+ DoFTools::extract_locally_relevant_dofs(dof_handler,
+ locally_relevant_dofs);
+
+ std::vector<types::global_dof_index> dofs_per_component(dim+1);
+ DoFTools::count_dofs_per_component(dof_handler, dofs_per_component);
+
+ // Discontinuous Galerkin methods are fantanistic methods in part because
+ // many of the limitations of traditional finite element methods no longer
+ // exist. Specifically, the need to use constraint matrices
+ // 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.
+ constraints.clear();
+ constraints.close();
+
+
+ // Just like step-40 we create a dynamic sparsity pattern
+ // and distribute it to the processors. Notice how we do not have to
+ // explictly mention that we are using a FESystem for system of
+ // variables instead of a FE_DGQ for a scalar variable
+ // or that we are using a discributed DoFHandler. All these specifics
+ // are taken care of under the hood by the deal.ii library.
+ // In order to build the sparsity
+ // pattern we use the DoFTools::make_flux_sparsity_pattern function
+ // since we using a DG method and need to take into account the DG
+ // fluxes in the sparsity pattern.
+ DynamicSparsityPattern dsp(dof_handler.n_dofs());
+ DoFTools::make_flux_sparsity_pattern(dof_handler,
+ dsp);
+
+ SparsityTools::distribute_sparsity_pattern(dsp,
+ dof_handler.n_locally_owned_dofs_per_processor(),
+ MPI_COMM_WORLD,
+ locally_relevant_dofs);
+
+ // Here is one area that I had to learn the hard way. The local
+ // discontinuous Galerkin method like the Mixed Method is written
+ // in mixed form and will lead to a block-structured matrix.
+ // In step-20 we see that we that we initialize the
+ // <code>system_martrix</code>
+ // such that we explicitly declare it to be block-structured.
+ // It turns out there are reasons to do this when you are going to be
+ // using a Schur complement method to solve the system of equations.
+ // While the LDG method will lead to a block-structured matrix,
+ // we do not have to explicitly declare our matrix to be one.
+ // I found that most of the distributed linear solvers did not
+ // accept block structured matrices and since I was using a
+ // distributed direct solver it was unnecessary to explicitly use a
+ // block structured matrix.
+ system_matrix.reinit(locally_owned_dofs,
+ locally_owned_dofs,
+ dsp,
+ MPI_COMM_WORLD);
+
+ // The final note that I will make in that this subroutine is that
+ // we initialize this processors solution and the
+ // right hand side vector the exact same was as we did in step-40.
+ // We should note that the <code>locally_relevant_solution</code> solution
+ // vector includes dofs that are locally relevant to our computations
+ // while the <code>system_rhs</code> right hand side vector will only
+ // include dofs that are locally owned by this processor.
+ locally_relevant_solution.reinit(locally_relevant_dofs,
+ MPI_COMM_WORLD);
+
+ system_rhs.reinit(locally_owned_dofs,
+ locally_relevant_dofs,
+ MPI_COMM_WORLD,
+ true);
+
+ const unsigned int n_vector_field = dim * dofs_per_component[0];
+ const unsigned int n_potential = dofs_per_component[1];
+
+ pcout << "Number of active cells : "
+ << triangulation.n_global_active_cells()
+ << std::endl
+ << "Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << " (" << n_vector_field << " + " << n_potential << ")"
+ << std::endl;
+}
+
+//@sect4{assemble_system}
+// This is the function that will assemble the global system matrix and
+// global right hand side vector for the LDG method. It starts out
+// like many of the deal.ii tutorial codes: declaring quadrature and
+// UpdateFlags objects, as well as vectors to that will hold the
+// dof indices for the cells we are working on in the global system.
+template <int dim>
+void
+LDGPoissonProblem<dim>::
+assemble_system()
+{
+ TimerOutput::Scope t(computing_timer, "assembly");
+
+ QGauss<dim> quadrature_formula(fe.degree+1);
+ QGauss<dim-1> face_quadrature_formula(fe.degree+1);
+
+ const UpdateFlags update_flags = update_values
+ | update_gradients
+ | update_quadrature_points
+ | update_JxW_values;
+
+ const UpdateFlags face_update_flags = update_values
+ | update_normal_vectors
+ | update_quadrature_points
+ | update_JxW_values;
+
+ 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);
+
+ // We first remark that we have the FEValues objects for
+ // the values of our cell basis functions as was done in most
+ // other examples. Now because we
+ // are using discontinuous Galerkin methods we also introduce a
+ // 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_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
+ // 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.
+ FEValues<dim> fe_values(fe, quadrature_formula, update_flags);
+
+ 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);
+
+ FESubfaceValues<dim> fe_subface_values(fe, face_quadrature_formula,
+ face_update_flags);
+
+ // Here are the local (dense) matrix and right hand side vector for
+ // the solid integrals as well as the integrals on the boundaries in the
+ // local discontinuous Galerkin method. These terms will be built for
+ // each local element in the mesh and then distributed to the global
+ // system matrix and right hand side vector.
+ FullMatrix<double> local_matrix(dofs_per_cell,dofs_per_cell);
+ Vector<double> local_vector(dofs_per_cell);
+
+ // The next four matrices are used to incorporate the flux integrals across
+ // interior faces of the mesh:
+ FullMatrix<double> vi_ui_matrix(dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> vi_ue_matrix(dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> ve_ui_matrix(dofs_per_cell, dofs_per_cell);
+ FullMatrix<double> ve_ue_matrix(dofs_per_cell, dofs_per_cell);
+ // As explained in the section on the LDG method we take our test
+ // function to be v and multiply it on the left side of our differential
+ // equation that is on u and peform integration by parts as explain in the
+ // introduction. Using this notation for test and solution function,
+ // the matrices below will then stand for:
+ //
+ // <code>vi_ui</code> - Taking the value of the test function from
+ // interior of this cell's face and the solution function
+ // from the interior of this cell.
+ //
+ // <code>vi_ue</code> - Taking the value of the test function from
+ // interior of this cell's face and the solution function
+ // from the exterior of this cell.
+ //
+ // <code>ve_ui</code> - Taking the value of the test function from
+ // exterior of this cell's face and the solution function
+ // from the interior of this cell.
+ //
+ // <code>ve_ue</code> - Taking the value of the test function from
+ // exterior of this cell's face and the solution function
+ // from the exterior of this cell.
+
+ // Now that we have gotten preliminary orders out of the way,
+ // we loop over all the cells
+ // and assemble the local system matrix and local right hand side vector
+ // using the DoFHandler::active_cell_iterator,
+ typename DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+
+ for(; cell!=endc; cell++)
+ {
+ // Now, since we are working in a distributed setting,
+ // we can only work on cells and write to dofs in the
+ // <code>system_matrix</code>
+ // and <code>rhs_vector</code>
+ // that corresponds to cells that are locally owned
+ // by this processor. We note that while we can only write to locally
+ // owned dofs, we will still use information from cells that are
+ // locally relevant. This is very much the same as in step-40.
+ if(cell->is_locally_owned())
+ {
+ // We now assemble the local contributions to the system matrix
+ // that includes the solid integrals in the LDG method as well as
+ // the right hand side vector. This involves resetting the local
+ // matrix and vector to contain all zeros, reinitializing the
+ // FEValues object for this cell and then building the
+ // <code>local_matrix</code> and <code>local_rhs</code> vector.
+ local_matrix = 0;
+ local_vector = 0;
+
+ fe_values.reinit(cell);
+ assemble_cell_terms(fe_values,
+ local_matrix,
+ local_vector);
+
+ // We remark that we need to get the local indices for the dofs to
+ // to this cell before we begin to compute the contributions
+ // from the numerical fluxes, i.e. the boundary conditions and
+ // interior fluxes.
+ cell->get_dof_indices(local_dof_indices);
+
+ // Now is where were start to loop over all the faces of the cell
+ // 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
+ // 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
+ // 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
+ // in the <code>local_matrix</code> as well as the boundary
+ // condition that ends up
+ // in the <code>local_vector</code>.
+ for(unsigned int face_no=0;
+ face_no< GeometryInfo<dim>::faces_per_cell;
+ face_no++)
+ {
+ typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+
+ if(face->at_boundary() )
+ {
+ fe_face_values.reinit(cell, face_no);
+
+ if(face->boundary_id() == Dirichlet)
+ {
+ // Notice here that in order to assemble the
+ // flux due to the penalty term for the the
+ // Dirichlet boundary condition we need the
+ // local cell diameter size and we can get
+ // that value for this specific cell with
+ // the following,
+ double h = cell->diameter();
+ assemble_Dirichlet_boundary_terms(fe_face_values,
+ local_matrix,
+ local_vector,
+ h);
+ }
+ else if(face->boundary_id() == Neumann)
+ {
+ assemble_Neumann_boundary_terms(fe_face_values,
+ local_matrix,
+ local_vector);
+ }
+ else
+ Assert(false, ExcNotImplemented() );
+ }
+ 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());
+
+ 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.
+ if(face->has_children())
+ {
+ // 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.
+ 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.
+ typename DoFHandler<dim>::cell_iterator neighbor_child =
+ cell->neighbor_child_on_subface(face_no, subface_no);
+
+ Assert(!neighbor_child->has_children(), ExcInternalError());
+
+ // 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
+ // <code>neighbor_child</code>'s
+ // FEFaceValues and the FESubfaceValues
+ // objects on the appropriate faces.
+ vi_ui_matrix = 0;
+ vi_ue_matrix = 0;
+ ve_ui_matrix = 0;
+ ve_ue_matrix = 0;
+
+ fe_subface_values.reinit(cell, face_no, subface_no);
+ 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.
+ 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.
+ assemble_flux_terms(fe_subface_values,
+ fe_neighbor_face_values,
+ vi_ui_matrix,
+ vi_ue_matrix,
+ ve_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.
+ 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
+ // <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);
+ }
+ }
+ 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,
+ // 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.
+ 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>.
+ const unsigned int neighbor_face_no =
+ cell->neighbor_of_neighbor(face_no);
+
+ vi_ui_matrix = 0;
+ vi_ue_matrix = 0;
+ ve_ui_matrix = 0;
+ ve_ue_matrix = 0;
+
+ fe_face_values.reinit(cell, 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
+ // <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>distribute_local_flux_to_global()</code>
+ assemble_flux_terms(fe_face_values,
+ fe_neighbor_face_values,
+ vi_ui_matrix,
+ vi_ue_matrix,
+ ve_ui_matrix,
+ ve_ue_matrix,
+ h);
+
+ 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);
+
+
+ }
+ }
+ }
+ }
+
+
+ // Now that have looped over all the faces for this
+ // cell and computed as well as disributed the local
+ // flux matrices to the <code>system_matrix</code>, we
+ // can finally distribute the cell's <code>local_matrix</code>
+ // and <code>local_vector</code> contribution to the
+ // global system matrix and global right hand side vector.
+ // We remark that we have to wait until this point
+ // to distribute the <code>local_matrix</code>
+ // and <code>system_rhs</code> to the global system.
+ // The reason being that in looping over the faces
+ // the faces on the boundary of the domain contribute
+ // to the <code>local_matrix</code>
+ // and <code>system_rhs</code>. We could distribute
+ // the local contributions for each component seperately,
+ // but writing to the distributed sparse matrix and vector
+ // is expensive and want to to minimize the number of times
+ // we do so.
+ constraints.distribute_local_to_global(local_matrix,
+ local_dof_indices,
+ system_matrix);
+
+ constraints.distribute_local_to_global(local_vector,
+ local_dof_indices,
+ system_rhs);
+
+ }
+ }
+
+ // We need to synchronize assembly of our global system
+ // matrix and global right hand side vector with all the other
+ // processors and
+ // use the compress() function to do this.
+ // This was discussed in detail in step-40.
+ system_matrix.compress(VectorOperation::add);
+ system_rhs.compress(VectorOperation::add);
+}
+
+
+// @sect4{assemble_cell_terms}
+// This function deals with constructing the local matrix due to
+// the solid integrals over each element and is very similar to the
+// the other examples in the deal.ii tutorials.
+template<int dim>
+void
+LDGPoissonProblem<dim>::
+assemble_cell_terms(
+ const FEValues<dim> &cell_fe,
+ FullMatrix<double> &cell_matrix,
+ Vector<double> &cell_vector)
+{
+ const unsigned int dofs_per_cell = cell_fe.dofs_per_cell;
+ const unsigned int n_q_points = cell_fe.n_quadrature_points;
+
+ const FEValuesExtractors::Vector VectorField(0);
+ const FEValuesExtractors::Scalar Potential(dim);
+
+ std::vector<double> rhs_values(n_q_points);
+
+ // We first get the value of the right hand side function
+ // evaluated at the quadrature points in the cell.
+ rhs_function.value_list(cell_fe.get_quadrature_points(),
+ rhs_values);
+
+ // Now, we loop over the quadrature points in the
+ // cell and then loop over the degrees of freedom and perform
+ // quadrature to approximate the integrals.
+ for(unsigned int q=0; q<n_q_points; q++)
+ {
+ 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);
+
+ for(unsigned int j=0; j<dofs_per_cell; j++)
+ {
+ const Tensor<1, dim> psi_j_field = cell_fe[VectorField].value(j,q);
+ const double psi_j_potential = cell_fe[Potential].value(j,q);
+
+ // This computation corresponds to assembling the local system
+ // matrix for the integral over an element,
+ //
+ // $\int_{\Omega_{e}} \left(\textbf{w} \cdot \textbf{q}
+ // - \nabla \cdot \textbf{w} u
+ // - \nabla w \cdot \textbf{q}
+ // \right) dx $
+ cell_matrix(i,j) += ( (psi_i_field * psi_j_field)
+ -
+ (div_psi_i_field * psi_j_potential)
+ -
+ (grad_psi_i_potential * psi_j_field)
+ ) * cell_fe.JxW(q);
+ }
+
+ // And this local right hand vector corresponds to the integral
+ // over the element cell,
+ //
+ // $ \int_{\Omega_{e}} w \, f(\textbf{x}) \, dx $
+ cell_vector(i) += psi_i_potential *
+ rhs_values[q] *
+ cell_fe.JxW(q);
+ }
+ }
+}
+
+// @sect4{assemble_Dirichlet_boundary_terms}
+// Here we have the function that builds the <code>local_matrix</code>
+// contribution
+// and local right hand side vector, <code>local_vector</code>
+// for the Dirichlet boundary condtions.
+template<int dim>
+void
+LDGPoissonProblem<dim>::
+assemble_Dirichlet_boundary_terms(
+ const FEFaceValues<dim> &face_fe,
+ FullMatrix<double> &local_matrix,
+ Vector<double> &local_vector,
+ const double & h)
+{
+ const unsigned int dofs_per_cell = face_fe.dofs_per_cell;
+ const unsigned int n_q_points = face_fe.n_quadrature_points;
+
+ const FEValuesExtractors::Vector VectorField(0);
+ const FEValuesExtractors::Scalar Potential(dim);
+
+ std::vector<double> Dirichlet_bc_values(n_q_points);
+
+ // In order to evaluate the flux on the Dirichlet boundary face we
+ // first get the value of the Dirichlet boundary function on the quadrature
+ // points of the face. Then we loop over all the quadrature points and
+ // degrees of freedom and approximate the integrals on the Dirichlet boundary
+ // element faces.
+ Dirichlet_bc_function.value_list(face_fe.get_quadrature_points(),
+ Dirichlet_bc_values);
+
+ for(unsigned int q=0; q<n_q_points; 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);
+
+ 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);
+
+ // We compute contribution for the flux $\widehat{q}$ on
+ // the Dirichlet boundary which enters our system matrix as,
+ //
+ // $ \int_{\text{face}} w \, ( \textbf{n} \cdot \textbf{q}
+ // + \sigma u) ds $
+ local_matrix(i,j) += psi_i_potential * (
+ face_fe.normal_vector(q) *
+ psi_j_field
+ +
+ (penalty/h) *
+ psi_j_potential) *
+ face_fe.JxW(q);
+
+ }
+
+ // We also compute the contribution for the flux for $\widehat{u}$
+ // on the Dirichlet boundary which is the Dirichlet boundary
+ // condition function and enters the right hand side vector as
+ //
+ // $\int_{\text{face}} (-\textbf{w} \cdot \textbf{n}
+ // + \sigma w) \, u_{D} ds $
+ local_vector(i) += (-1.0 * psi_i_field *
+ face_fe.normal_vector(q)
+ +
+ (penalty/h) *
+ psi_i_potential) *
+ Dirichlet_bc_values[q] *
+ face_fe.JxW(q);
+ }
+ }
+}
+
+// @sect4{assemble_Neumann_boundary_terms}
+// Here we have the function that builds the <code>local_matrix</code>
+// and <code>local_vector</code> for the Neumann boundary condtions.
+template<int dim>
+void
+LDGPoissonProblem<dim>::
+assemble_Neumann_boundary_terms(
+ const FEFaceValues<dim> &face_fe,
+ FullMatrix<double> &local_matrix,
+ Vector<double> &local_vector)
+{
+ const unsigned int dofs_per_cell = face_fe.dofs_per_cell;
+ const unsigned int n_q_points = face_fe.n_quadrature_points;
+
+ const FEValuesExtractors::Vector VectorField(0);
+ const FEValuesExtractors::Scalar Potential(dim);
+
+ // In order to get evaluate the flux on the Neumann boundary face we
+ // first get the value of the Neumann boundary function on the quadrature
+ // points of the face. Then we loop over all the quadrature points and
+ // degrees of freedom and approximate the integrals on the Neumann boundary
+ // element faces.
+ std::vector<double > Neumann_bc_values(n_q_points);
+
+ for(unsigned int q=0; q<n_q_points; 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);
+
+ for(unsigned int j=0; j<dofs_per_cell; j++)
+ {
+
+ 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,
+ //
+ // $\int_{\text{face}} \textbf{w} \cdot \textbf{n} \, u \, ds $
+ local_matrix(i,j) += psi_i_field *
+ face_fe.normal_vector(q) *
+ psi_j_potential *
+ face_fe.JxW(q);
+
+ }
+
+ // We also compute the contribution for the flux for
+ // $\widehat{q}$ on the Neumann bounary which is the
+ // Neumann boundary condition and enters the right
+ // hand side vector as
+ //
+ // $\int_{\text{face}} -w \, g_{N} \, ds$
+ local_vector(i) += -psi_i_potential *
+ Neumann_bc_values[q] *
+ face_fe.JxW(q);
+ }
+ }
+}
+
+// @sect4{assemble_flux_terms}
+// Now we finally get to the function which builds the interior fluxes.
+// This is a rather long function
+// and we will describe what is going on in detail.
+template<int dim>
+void
+LDGPoissonProblem<dim>::
+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)
+{
+ const unsigned int n_face_points = fe_face_values.n_quadrature_points;
+ const unsigned int dofs_this_cell = fe_face_values.dofs_per_cell;
+ const unsigned int dofs_neighbor_cell = fe_neighbor_face_values.dofs_per_cell;
+
+ const FEValuesExtractors::Vector VectorField(0);
+ const FEValuesExtractors::Scalar Potential(dim);
+
+ // The first thing we do is after the boilerplate is define
+ // the unit vector $\boldsymbol \beta$ that is used in defining
+ // the LDG/ALternating fluxes.
+ Point<dim> beta;
+ for(int i=0; i<dim; i++)
+ beta(i) = 1.0;
+ beta /= sqrt(beta.square() );
+
+ // Now we loop over all the quadrature points on the element face
+ // and loop over all the degrees of freedom and approximate
+ // the following flux integrals.
+ for(unsigned int q=0; q<n_face_points; q++)
+ {
+ for(unsigned int i=0; i<dofs_this_cell; i++)
+ {
+ const Tensor<1,dim> psi_i_field_minus =
+ fe_face_values[VectorField].value(i,q);
+ const double psi_i_potential_minus =
+ fe_face_values[Potential].value(i,q);
+
+ for(unsigned int j=0; j<dofs_this_cell; j++)
+ {
+ const Tensor<1,dim> psi_j_field_minus =
+ fe_face_values[VectorField].value(j,q);
+ const double psi_j_potential_minus =
+ fe_face_values[Potential].value(j,q);
+
+ // We compute the flux matrix where the test function's
+ // as well as the solution function's values are taken from
+ // the interior as,
+ //
+ // $\int_{\text{face}}
+ // \left( \frac{1}{2} \, n^{-} \cdot ( \textbf{w}^{-} u^{-}
+ // + w^{-} \textbf{q}^{-})
+ // + \boldsymbol \beta \cdot \textbf{w}^{-} u^{-}
+ // - w^{-} \boldsymbol \beta \cdot \textbf{q}^{-}
+ // + \sigma w^{-} \, u^{-} \right) ds$
+ vi_ui_matrix(i,j) += (0.5 * (
+ psi_i_field_minus *
+ fe_face_values.normal_vector(q) *
+ psi_j_potential_minus
+ +
+ psi_i_potential_minus *
+ fe_face_values.normal_vector(q) *
+ psi_j_field_minus )
+ +
+ beta *
+ psi_i_field_minus *
+ psi_j_potential_minus
+ -
+ beta *
+ psi_i_potential_minus *
+ psi_j_field_minus
+ +
+ (penalty/h) *
+ psi_i_potential_minus *
+ psi_j_potential_minus
+ ) *
+ fe_face_values.JxW(q);
+ }
+
+ for(unsigned int j=0; j<dofs_neighbor_cell; j++)
+ {
+ const Tensor<1,dim> psi_j_field_plus =
+ fe_neighbor_face_values[VectorField].value(j,q);
+ const double psi_j_potential_plus =
+ fe_neighbor_face_values[Potential].value(j,q);
+
+ // We compute the flux matrix where the test function is
+ // from the interior of this elements face and solution
+ // function is taken from the exterior. This corresponds
+ // to the computation,
+ //
+ // $\int_{\text{face}}
+ // \left( \frac{1}{2} \, n^{-} \cdot ( \textbf{w}^{-} u^{+}
+ // + w^{-} \textbf{q}^{+})
+ // - \boldsymbol \beta \cdot \textbf{w}^{-} u^{+}
+ // + w^{-} \boldsymbol \beta \cdot \textbf{q}^{+}
+ // - \sigma w^{-} \, u^{+} \right) ds $
+ vi_ue_matrix(i,j) += ( 0.5 * (
+ psi_i_field_minus *
+ fe_face_values.normal_vector(q) *
+ psi_j_potential_plus
+ +
+ psi_i_potential_minus *
+ fe_face_values.normal_vector(q) *
+ psi_j_field_plus )
+ -
+ beta *
+ psi_i_field_minus *
+ psi_j_potential_plus
+ +
+ beta *
+ psi_i_potential_minus *
+ psi_j_field_plus
+ -
+ (penalty/h) *
+ psi_i_potential_minus *
+ psi_j_potential_plus
+ ) *
+ fe_face_values.JxW(q);
+ }
+ }
+
+ for(unsigned int i=0; i<dofs_neighbor_cell; i++)
+ {
+ const Tensor<1,dim> psi_i_field_plus =
+ fe_neighbor_face_values[VectorField].value(i,q);
+ const double psi_i_potential_plus =
+ fe_neighbor_face_values[Potential].value(i,q);
+
+ for(unsigned int j=0; j<dofs_this_cell; j++)
+ {
+ const Tensor<1,dim> psi_j_field_minus =
+ fe_face_values[VectorField].value(j,q);
+ const double psi_j_potential_minus =
+ fe_face_values[Potential].value(j,q);
+
+
+ // We compute the flux matrix where the test function is
+ // from the exterior of this elements face and solution
+ // function is taken from the interior. This corresponds
+ // to the computation,
+ //
+ // $ \int_{\text{face}}
+ // \left( -\frac{1}{2}\, n^{-} \cdot (\textbf{w}^{+} u^{-}
+ // + w^{+} \textbf{q}^{-} )
+ // - \boldsymbol \beta \cdot \textbf{w}^{+} u^{-}
+ // + w^{+} \boldsymbol \beta \cdot \textbf{q}^{-}
+ // - \sigma w^{+} u^{-} \right) ds $
+ ve_ui_matrix(i,j) += (-0.5 * (
+ psi_i_field_plus *
+ fe_face_values.normal_vector(q) *
+ psi_j_potential_minus
+ +
+ psi_i_potential_plus *
+ fe_face_values.normal_vector(q) *
+ psi_j_field_minus)
+ -
+ beta *
+ psi_i_field_plus *
+ psi_j_potential_minus
+ +
+ beta *
+ psi_i_potential_plus *
+ psi_j_field_minus
+ -
+ (penalty/h) *
+ psi_i_potential_plus *
+ psi_j_potential_minus
+ ) *
+ fe_face_values.JxW(q);
+ }
+
+ for(unsigned int j=0; j<dofs_neighbor_cell; j++)
+ {
+ const Tensor<1,dim> psi_j_field_plus =
+ fe_neighbor_face_values[VectorField].value(j,q);
+ const double psi_j_potential_plus =
+ fe_neighbor_face_values[Potential].value(j,q);
+
+ // And lastly we compute the flux matrix where the test
+ // function and solution function are taken from the exterior
+ // cell to this face. This corresponds to the computation,
+ //
+ // $\int_{\text{face}}
+ // \left( -\frac{1}{2}\, n^{-} \cdot ( \textbf{w}^{+} u^{+}
+ // + w^{+} \textbf{q}^{+} )
+ // + \boldsymbol \beta \cdot \textbf{w}^{+} u^{+}
+ // - w^{+} \boldsymbol \beta \cdot \textbf{q}^{+}
+ // + \sigma w^{+} u^{+} \right) ds $
+ ve_ue_matrix(i,j) += (-0.5 * (
+ psi_i_field_plus *
+ fe_face_values.normal_vector(q) *
+ psi_j_potential_plus
+ +
+ psi_i_potential_plus *
+ fe_face_values.normal_vector(q) *
+ psi_j_field_plus )
+ +
+ beta *
+ psi_i_field_plus *
+ psi_j_potential_plus
+ -
+ beta *
+ psi_i_potential_plus *
+ psi_j_field_plus
+ +
+ (penalty/h) *
+ psi_i_potential_plus *
+ psi_j_potential_plus
+ ) *
+ fe_face_values.JxW(q);
+ }
+
+ }
+ }
+}
+
+// @sect4{distribute_local_flux_to_global}
+// In this function we use the ConstraintMatrix to distribute
+// the local flux matrices to the global system matrix.
+// Since I have to do this twice in assembling the
+// system matrix, I made function to do it rather than have
+// repeated code.
+// We remark that the reader take special note of
+// the which matrices we are distributing and the order
+// in which we pass the dof indices vectors. In distributing
+// the first matrix, i.e. <code>vi_ui_matrix</code>, we are
+// taking the test function and solution function values from
+// the interior of this cell and therefore only need the
+// <code>local_dof_indices</code> since it contains the dof
+// indices to this cell. When we distribute the second matrix,
+// <code>vi_ue_matrix</code>, the test function is taken
+// form the inteior of
+// this cell while the solution function is taken from the
+// exterior, i.e. the neighbor cell. Notice that the order
+// degrees of freedom index vectors matrch this pattern: first
+// the <code>local_dof_indices</code> which is local to
+// this cell and then
+// the <code>local_neighbor_dof_indices</code> which is
+// local to the neighbor's
+// cell. The order in which we pass the dof indices for the
+// matrices is paramount to constructing our global system
+// matrix properly. The ordering of the last two matrices
+// follow the same logic as the first two we discussed.
+template<int dim>
+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)
+{
+ constraints.distribute_local_to_global(vi_ui_matrix,
+ local_dof_indices,
+ system_matrix);
+
+ constraints.distribute_local_to_global(vi_ue_matrix,
+ local_dof_indices,
+ local_neighbor_dof_indices,
+ system_matrix);
+
+ constraints.distribute_local_to_global(ve_ui_matrix,
+ local_neighbor_dof_indices,
+ local_dof_indices,
+ system_matrix);
+
+ constraints.distribute_local_to_global(ve_ue_matrix,
+ local_neighbor_dof_indices,
+ system_matrix);
+}
+
+
+
+// @sect4{solve}
+// As mentioned earlier I used a direct solver to solve
+// the linear system of equations resulting from the LDG
+// method applied to the Poisson equation. One could also
+// use a iterative sovler, however, we then need to use
+// a preconditoner and that was something I did not wanted
+// to get into. The uses of a direct sovler here is
+// somewhat of a limitation. The built-in distributed
+// direct solver in Trilinos reduces everything to one
+// processor, solves the system and then distributes
+// everything back out to the other processors. However,
+// by linking to more advanced direct sovlers through
+// Trilinos one can accomplish fully distributed computations
+// and not much about the following function calls will
+// change.
+template<int dim>
+void
+LDGPoissonProblem<dim>::
+solve()
+{
+ TimerOutput::Scope t(computing_timer, "solve");
+
+ // As in step-40 in order to perform a linear solve
+ // we need solution vector where there is no overlap across
+ // the processors and we create this by instantiating
+ // <code>completely_distributed_solution</code> solution
+ // vector using
+ // the copy constructor on the global system right hand
+ // side vector which itself is completely distributed vector.
+ TrilinosWrappers::MPI::Vector
+ completely_distributed_solution(system_rhs);
+
+ // Now we can preform the solve on the completeley distributed
+ // right hand side vector, system matrix and the completely
+ // distributed solution.
+ solver.solve(system_matrix,
+ completely_distributed_solution,
+ system_rhs);
+
+ // We now distribute the constraints of our system onto the
+ // completely solution vector, but in our case with the LDG
+ // method there are none.
+ constraints.distribute(completely_distributed_solution);
+
+ // Lastly we copy the completely distributed solution vector,
+ // <code>completely_distributed_solution</code>,
+ // to solution vector which has some overlap between
+ // processors, <code>locally_relevant_solution</code>.
+ // We need the overlapped portions of our solution
+ // in order to be able to do computations using the solution
+ // later in the code or in post processing.
+ locally_relevant_solution = completely_distributed_solution;
+}
+
+// @sect4{output_results}
+// This function deals with the writing of the reuslts in parallel
+// to disk. It is almost exactly the same as
+// in step-40 and we wont go into it. It is noteworthy
+// that in step-40 the output is only the scalar solution,
+// while in our situation, we are outputing both the scalar
+// solution as well as the vector field solution. The only
+// difference between this function and the one in step-40
+// is in the <code>solution_names</code> vector where we have to add
+// the gradient dimensions. Everything else is taken care
+// of by the deal.ii library!
+template<int dim>
+void
+LDGPoissonProblem<dim>::
+output_results() const
+{
+ std::vector<std::string> solution_names;
+ switch(dim)
+ {
+ case 1:
+ solution_names.push_back("u");
+ solution_names.push_back("du/dx");
+ break;
+
+ case 2:
+ solution_names.push_back("grad(u)_x");
+ solution_names.push_back("grad(u)_y");
+ solution_names.push_back("u");
+ break;
+
+ case 3:
+ solution_names.push_back("grad(u)_x");
+ solution_names.push_back("grad(u)_y");
+ solution_names.push_back("grad(u)_z");
+ solution_names.push_back("u");
+ break;
+
+ default:
+ Assert(false, ExcNotImplemented() );
+ }
+
+ DataOut<dim> data_out;
+ data_out.attach_dof_handler(dof_handler);
+ data_out.add_data_vector(locally_relevant_solution,
+ solution_names);
+
+ Vector<float> subdomain(triangulation.n_active_cells());
+
+ for(unsigned int i=0; i<subdomain.size(); i++)
+ subdomain(i) = triangulation.locally_owned_subdomain();
+
+ data_out.add_data_vector(subdomain,"subdomain");
+
+ data_out.build_patches();
+
+ const std::string filename = ("solution." +
+ Utilities::int_to_string(
+ triangulation.locally_owned_subdomain(),4));
+
+ std::ofstream output((filename + ".vtu").c_str());
+ data_out.write_vtu(output);
+
+ if(Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0 )
+ {
+ std::vector<std::string> filenames;
+ for(unsigned int i=0;
+ i < Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+ i++)
+ {
+ filenames.push_back("solution." +
+ Utilities::int_to_string(i,4) +
+ ".vtu");
+ }
+ std::ofstream master_output("solution.pvtu");
+ data_out.write_pvtu_record(master_output, filenames);
+ }
+}
+
+
+// @sect4{run}
+// The only public function of this class is pretty much exactly
+// the same as all the other deal.ii examples except I setting
+// the constant in the DG penalty ($\tilde{\sigma}}$) to be 1.
+template<int dim>
+void
+LDGPoissonProblem<dim>::
+run()
+{
+ penalty = 1.0;
+ make_grid();
+ make_dofs();
+ assemble_system();
+ solve();
+ output_results();
+}
+
+
+// @sect3{main}
+// Here it the main class of our program, since it is nearly exactly
+// the same as step-40 and many of the other examples I won't
+// elaborate on it.
+int main(int argc, char *argv[])
+{
+
+ try {
+ using namespace dealii;
+
+ deallog.depth_console(0);
+
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv,
+ numbers::invalid_unsigned_int);
+
+ unsigned int degree = 1;
+ unsigned int n_refine = 6;
+ LDGPoissonProblem<2> Poisson(degree, n_refine);
+ Poisson.run();
+
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+ return 0;
+}
projects / code-gallery.git / commitdiff