vector_fe_ex7.C File Reference

Go to the source code of this file.

Typedefs
typedef MatrixXcd	EigenMatrix
typedef VectorXcd	EigenVector

Functions
void	fe_assembly (EquationSystems &es, bool global_solve)
void	assemble_divgrad (EquationSystems &es, const std::string &system_name)
int	main (int argc, char **argv)
void	assemble_divgrad (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
int	main ()

Typedef Documentation

EigenMatrix

typedef MatrixXd EigenMatrix

Definition at line 149 of file vector_fe_ex7.C.

EigenVector

typedef VectorXd EigenVector

Definition at line 150 of file vector_fe_ex7.C.

Function Documentation

assemble_divgrad() [1/2]

void assemble_divgrad	(	EquationSystems &	es,
		const std::string &	libmesh_dbg_varsystem_name
	)

Definition at line 726 of file vector_fe_ex7.C.

{
  libmesh_assert_equal_to(system_name, "Lambda");
  fe_assembly(es, /*global_solve=*/true);
}

References fe_assembly().

assemble_divgrad() [2/2]

void assemble_divgrad	(	EquationSystems &	es,
		const std::string &	system_name
	)

Referenced by main().

fe_assembly()

void fe_assembly	(	EquationSystems &	es,
		bool	global_solve
	)

The scalar solution at the quadrature points

Data structures for computing the enriched scalar solution

Definition at line 326 of file vector_fe_ex7.C.

{
  const MeshBase & mesh = es.get_mesh();
  const unsigned int dim = mesh.mesh_dimension();
  // Are our elements simplicial?
  const bool simplicial = es.parameters.get<bool>("simplicial");
 
  // The div-grad, e.g. vector-scalar system
  auto & system = es.get_system<System>("DivGrad");
  // Our implicit Lagrange multiplier system
  auto & lambda_system = es.get_system<LinearImplicitSystem>("Lambda");
 
  const auto & dof_map = system.get_dof_map();
  const auto & lambda_dof_map = lambda_system.get_dof_map();
 
  const FEType vector_fe_type = dof_map.variable_type(system.variable_number("u"));
  const FEType scalar_fe_type = dof_map.variable_type(system.variable_number("p"));
  const FEType enriched_scalar_fe_type =
      dof_map.variable_type(system.variable_number("p_enriched"));
  const FEType lambda_fe_type =
      lambda_dof_map.variable_type(lambda_system.variable_number("lambda"));
 
  // Volumetric FE objects
  std::unique_ptr<FEVectorBase> vector_fe(FEVectorBase::build(dim, vector_fe_type));
  std::unique_ptr<FEBase> scalar_fe(FEBase::build(dim, scalar_fe_type));
  std::unique_ptr<FEBase> enriched_scalar_fe(FEBase::build(dim, enriched_scalar_fe_type));
 
  // Volumetric quadrature rule
  QGauss qrule(dim, FIFTH);
 
  // Attach quadrature rules for the FE objects that we will reinit within the element "volume"
  vector_fe->attach_quadrature_rule(&qrule);
  scalar_fe->attach_quadrature_rule(&qrule);
  enriched_scalar_fe->attach_quadrature_rule(&qrule);
 
  // Declare finite element objects for boundary integration
  std::unique_ptr<FEVectorBase> vector_fe_face(FEVectorBase::build(dim, vector_fe_type));
  std::unique_ptr<FEBase> enriched_scalar_fe_face(FEBase::build(dim, enriched_scalar_fe_type));
  std::unique_ptr<FEBase> lambda_fe_face(FEBase::build(dim, lambda_fe_type));
 
  // Boundary integration requires one quadrature rule with dimensionality one
  // less than the dimensionality of the element.
  QGauss qface(dim - 1, FIFTH);
 
  // Attach quadrature rules for the FE objects that we will reinit on the element faces
  vector_fe_face->attach_quadrature_rule(&qface);
  enriched_scalar_fe_face->attach_quadrature_rule(&qface);
  lambda_fe_face->attach_quadrature_rule(&qface);
 
  // pre-request our required volumetric data
  const auto & JxW = vector_fe->get_JxW();
  const auto & q_point = vector_fe->get_xyz();
  const auto & vector_phi = vector_fe->get_phi();
  const auto & scalar_phi = scalar_fe->get_phi();
  const auto & enriched_scalar_phi = enriched_scalar_fe->get_phi();
  const auto & div_vector_phi = vector_fe->get_div_phi();
 
  // pre-request our required element face data
  const auto & vector_phi_face = vector_fe_face->get_phi();
  const auto & enriched_scalar_phi_face = enriched_scalar_fe_face->get_phi();
  const auto & lambda_phi_face = lambda_fe_face->get_phi();
  const auto & JxW_face = vector_fe_face->get_JxW();
  const auto & qface_point = vector_fe_face->get_xyz();
  const auto & normals = vector_fe_face->get_normals();
 
  //
  // We follow the notation of Cockburn in "A Charactrization of Hybridized
  // Mixed methods for Second Order Elliptic problems"for the element
  // matrices/vectors. We will need "Eigen" versions of many of the
  // matrices/vectors because the libMesh DenseMatrix doesn't have an inverse
  // API
  //
 
  // LM matrix and RHS after eliminating vector and scalar dofs
  DenseMatrix<Number> E_libmesh;
  DenseVector<Number> H_libmesh;
  EigenMatrix E;
  EigenMatrix H;
  // Auxiliary matrices and RHS. A is vector-vector. B is vector-scalar. C is
  // vector-LM. Sinv is the inverse of the Schur complement S which is given by
  // S = Bt * A^{-1} * B. G is the RHS of the vector equation (e.g. the
  // Dirichlet condition). F is the RHS of the scalar equation (resulting from
  // the body force or MMS force in this case). Hg is the post-elimination
  // version of G, and Hf is the post-elimination version of F
  EigenMatrix A, Ainv, B, Bt, C, Ct, Sinv;
  EigenVector G, F, Hg, Hf;
  // element matrix for boundary LM dofs. We have to have this because our
  // SIDE_HIERARCHIC variables live on *all* element faces, not just interior
  // ones.
  EigenMatrix L;
 
  // Lambda eigen vector for constructing vector and scalar solutions
  EigenVector Lambda;
  // The lambda solution at the quadrature points
  std::vector<Number> lambda_qps;
  std::vector<Number> scalar_qps;
 
  DenseMatrix<Number> K_enriched_scalar;
  DenseVector<Number> F_enriched_scalar, U_enriched_scalar;
 
  // Containers for dof indices
  std::vector<dof_id_type> vector_dof_indices;
  std::vector<dof_id_type> scalar_dof_indices;
  std::vector<dof_id_type> enriched_scalar_dof_indices;
  std::vector<dof_id_type> lambda_dof_indices;
  std::vector<Number> lambda_solution_std_vec;
 
  // The global system matrix
  SparseMatrix<Number> & matrix = lambda_system.get_system_matrix();
 
  for (const auto & elem : mesh.active_local_element_ptr_range())
  {
    // Retrive our dof indices for all fields
    dof_map.dof_indices(elem, vector_dof_indices, system.variable_number("u"));
    dof_map.dof_indices(elem, scalar_dof_indices, system.variable_number("p"));
    lambda_dof_map.dof_indices(elem, lambda_dof_indices, lambda_system.variable_number("lambda"));
 
    const auto vector_n_dofs = vector_dof_indices.size();
    const auto scalar_n_dofs = scalar_dof_indices.size();
    const auto lambda_n_dofs = lambda_dof_indices.size();
 
    // Reinit our volume FE objects
    vector_fe->reinit(elem);
    scalar_fe->reinit(elem);
 
    libmesh_assert_equal_to(vector_n_dofs, vector_phi.size());
    libmesh_assert_equal_to(scalar_n_dofs, scalar_phi.size());
 
    // prepare our matrix/vector data structures for the vector equation
    A.setZero(vector_n_dofs, vector_n_dofs);
    B.setZero(vector_n_dofs, scalar_n_dofs);
    C.setZero(vector_n_dofs, lambda_n_dofs);
    G.setZero(vector_n_dofs);
    // and for the scalar equation
    Bt.setZero(scalar_n_dofs, vector_n_dofs);
    F.setZero(scalar_n_dofs);
    // and for the LM equation
    Ct.setZero(lambda_n_dofs, vector_n_dofs);
    L.setZero(lambda_n_dofs, lambda_n_dofs);
 
    for (const auto qp : make_range(qrule.n_points()))
    {
      // Vector equation dependence on vector dofs
      for (const auto i : make_range(vector_n_dofs))
        for (const auto j : make_range(vector_n_dofs))
          A(i, j) += JxW[qp] * (vector_phi[i][qp] * vector_phi[j][qp]);
 
      // Vector equation dependence on scalar dofs
      for (const auto i : make_range(vector_n_dofs))
        for (const auto j : make_range(scalar_n_dofs))
          B(i, j) -= JxW[qp] * (div_vector_phi[i][qp] * scalar_phi[j][qp]);
 
      // Scalar equation dependence on vector dofs
      for (const auto i : make_range(scalar_n_dofs))
        for (const auto j : make_range(vector_n_dofs))
          Bt(i, j) -= JxW[qp] * (scalar_phi[i][qp] * div_vector_phi[j][qp]);
 
      // This is the end of the matrix summation loop
      // Now we build the element right-hand-side contribution.
      // This involves a single loop in which we integrate the "forcing
      // function" in the PDE against the scalar test functions (k).
      {
        const Real x = q_point[qp](0);
        const Real y = q_point[qp](1);
        const Real z = q_point[qp](2);
 
        // "f" is the forcing function for the Poisson equation, which is
        // just the divergence of the exact solution for the vector field.
        // This is the well-known "method of manufactured solutions".
        Real f = 0;
        if (dim == 2)
          f = DivGradExactSolution().forcing(x, y);
        else if (dim == 3)
          f = DivGradExactSolution().forcing(x, y, z);
 
        // Scalar equation RHS
        for (const auto i : make_range(scalar_n_dofs))
          F(i) -= JxW[qp] * scalar_phi[i][qp] * f;
      }
    }
 
    {
      for (auto side : elem->side_index_range())
      {
        // Reinit our face FE objects
        vector_fe_face->reinit(elem, side);
        libmesh_assert_equal_to(vector_n_dofs, vector_phi_face.size());
        lambda_fe_face->reinit(elem, side);
 
        if (elem->neighbor_ptr(side) == nullptr)
          for (const auto qp : make_range(qface.n_points()))
          {
            const Real xf = qface_point[qp](0);
            const Real yf = qface_point[qp](1);
            const Real zf = qface_point[qp](2);
 
            // The boundary value for scalar field.
            Real scalar_value = 0;
            if (dim == 2)
              scalar_value = DivGradExactSolution().scalar(xf, yf);
            else if (dim == 3)
              scalar_value = DivGradExactSolution().scalar(xf, yf, zf);
 
            // External boundary -> Dirichlet faces -> Vector equaton RHS
            for (const auto i : make_range(vector_n_dofs))
              G(i) -= JxW_face[qp] * (vector_phi_face[i][qp] * normals[qp]) * scalar_value;
 
            // Need to do something with the external boundary LM dofs to prevent the matrix from
            // being singular
            for (const auto i : make_range(lambda_n_dofs))
              for (const auto j : make_range(lambda_n_dofs))
                L(i, j) += JxW_face[qp] * lambda_phi_face[i][qp] * lambda_phi_face[j][qp];
          }
        else
          for (const auto qp : make_range(qface.n_points()))
          {
            // Vector equation dependence on LM dofs
            for (const auto i : make_range(vector_n_dofs))
              for (const auto j : make_range(lambda_n_dofs))
                C(i, j) +=
                    JxW_face[qp] * (vector_phi_face[i][qp] * normals[qp]) * lambda_phi_face[j][qp];
 
            // LM equation dependence on vector dofs
            for (const auto i : make_range(lambda_n_dofs))
              for (const auto j : make_range(vector_n_dofs))
                Ct(i, j) +=
                    JxW_face[qp] * lambda_phi_face[i][qp] * (vector_phi_face[j][qp] * normals[qp]);
          }
      }
    }
 
    Ainv = A.inverse();
    // Compute the Schur complement inverse
    Sinv = (Bt * Ainv * B).inverse();
    // These equations can be derived at by eliminating the u and p dofs from the system
    // (A  B  C)(u)        (G)
    // (Bt 0  0)(p)      = (F)
    // (Ct 0  0)(lambda)   (0)
    E = Ct * Ainv * (A - B * Sinv * Bt) * Ainv * C;
    E = E + L;
    Hg = Ct * Ainv * (A - B * Sinv * Bt) * Ainv * G;
    Hf = Ct * Ainv * B * Sinv * F;
    H = Hg + Hf;
 
    // Build our libMesh data structures from the Eigen ones
    E_libmesh.resize(E.rows(), E.cols());
    H_libmesh.resize(H.size());
    libmesh_assert((E.rows() == E.cols()) && (E.rows() == H.size()));
    for (const auto i : make_range(E.rows()))
    {
      H_libmesh(i) = H(i);
      for (const auto j : make_range(E.cols()))
        E_libmesh(i, j) = E(i, j);
    }
 
    if (global_solve)
    {
      // We were performing our finite element assembly for the implicit solve step of our
      // example. Add our local element vectors/matrices into the global system
      dof_map.constrain_element_matrix_and_vector(E_libmesh, H_libmesh, lambda_dof_indices);
      matrix.add_matrix(E_libmesh, lambda_dof_indices);
      lambda_system.rhs->add_vector(H_libmesh, lambda_dof_indices);
    }
    else
    {
      //
      // We are doing our finite element assembly for the second time. We now know the Lagrange
      // multiplier solution. With that and the local element matrices and vectors we can compute
      // the vector and scalar solutions
      //
 
      Lambda.resize(lambda_n_dofs);
      lambda_system.current_local_solution->get(lambda_dof_indices, lambda_solution_std_vec);
      for (const auto i : make_range(lambda_n_dofs))
        Lambda(i) = lambda_solution_std_vec[i];
      const auto scalar_soln = Sinv * Bt * Ainv * G - Sinv * F - Sinv * Bt * Ainv * C * Lambda;
      const auto vector_soln = Ainv * (G - B * scalar_soln - C * Lambda);
      for (const auto i : make_range(vector_n_dofs))
        system.solution->set(vector_dof_indices[i], vector_soln(i));
      for (const auto i : make_range(scalar_n_dofs))
        system.solution->set(scalar_dof_indices[i], scalar_soln(i));
 
#if !defined(LIBMESH_HAVE_PETSC) || !defined(LIBMESH_USE_REAL_NUMBERS)
      if (!simplicial)
        // We don't support SVD solves in this configuration
        continue;
#endif
 
      //
      // Now solve for the enriched scalar solution using our Lagrange multiplier solution and, for
      // non-simplexes, the lower-order scalar solution. Note that the Lagrange multiplier
      // represents the trace of p so it is a logical choice to leverage in this postprocessing
      // stage!
      //
 
      dof_map.dof_indices(elem, enriched_scalar_dof_indices, system.variable_number("p_enriched"));
      const auto enriched_scalar_n_dofs = enriched_scalar_dof_indices.size();
      const auto m = simplicial ? lambda_n_dofs : lambda_n_dofs + 1;
      const auto n = enriched_scalar_n_dofs;
 
      K_enriched_scalar.resize(m, n);
      F_enriched_scalar.resize(m);
      U_enriched_scalar.resize(n);
 
      // L2 projection of the enriched scalar into the LM space
      for (const auto side : elem->side_index_range())
      {
        vector_fe_face->reinit(elem, side); // for JxW_face and qface_point
        enriched_scalar_fe_face->reinit(elem, side);
        lambda_fe_face->reinit(elem, side);
        if (elem->neighbor_ptr(side) == nullptr)
          for (const auto qp : make_range(qface.n_points()))
          {
            const Real xf = qface_point[qp](0);
            const Real yf = qface_point[qp](1);
            const Real zf = qface_point[qp](2);
 
            // The boundary value for scalar field.
            Real scalar_value = 0;
            if (dim == 2)
              scalar_value = DivGradExactSolution().scalar(xf, yf);
            else if (dim == 3)
              scalar_value = DivGradExactSolution().scalar(xf, yf, zf);
 
            for (const auto i : make_range(lambda_n_dofs))
            {
              F_enriched_scalar(i) += JxW_face[qp] * lambda_phi_face[i][qp] * scalar_value;
              for (const auto j : make_range(enriched_scalar_n_dofs))
                K_enriched_scalar(i, j) +=
                    JxW_face[qp] * lambda_phi_face[i][qp] * enriched_scalar_phi_face[j][qp];
            }
          }
        else
        {
          // compute local face lambda solution
          lambda_qps.resize(qface.n_points());
          for (auto & lambda_qp : lambda_qps)
            lambda_qp = 0;
          for (const auto qp : make_range(qface.n_points()))
            for (const auto i : index_range(lambda_solution_std_vec))
              lambda_qps[qp] += lambda_solution_std_vec[i] * lambda_phi_face[i][qp];
 
          for (const auto qp : make_range(qface.n_points()))
            for (const auto i : make_range(lambda_n_dofs))
            {
              F_enriched_scalar(i) += JxW_face[qp] * lambda_phi_face[i][qp] * lambda_qps[qp];
              for (const auto j : make_range(enriched_scalar_n_dofs))
                K_enriched_scalar(i, j) +=
                    JxW_face[qp] * lambda_phi_face[i][qp] * enriched_scalar_phi_face[j][qp];
            }
        }
      }
 
      if (simplicial)
        K_enriched_scalar.lu_solve(F_enriched_scalar, U_enriched_scalar);
      else
      {
        //
        // For tensor product elements, the system of equations coming from the L2 projection into
        // the Lagrange multiplier space is singular. To make the system nonsingular, we add an L2
        // projection into the low-order scalar solution space.
        //
 
        enriched_scalar_fe->reinit(elem);
        libmesh_assert(scalar_n_dofs == 1);
        libmesh_assert(scalar_soln.size() == 1);
        libmesh_assert(scalar_phi.size() == 1);
 
        scalar_qps.resize(qrule.n_points());
        for (auto & scalar_qp : scalar_qps)
          scalar_qp = 0;
        for (const auto qp : make_range(qrule.n_points()))
          scalar_qps[qp] += scalar_soln(0) * scalar_phi[0][qp];
        for (const auto qp : make_range(qrule.n_points()))
        {
          F_enriched_scalar(n) += JxW[qp] * scalar_phi[0][qp] * scalar_qps[qp];
          for (const auto j : make_range(n))
            K_enriched_scalar(n, j) += JxW[qp] * scalar_phi[0][qp] * enriched_scalar_phi[j][qp];
        }
        // We have more rows than columns so an LU solve isn't an option
        K_enriched_scalar.svd_solve(F_enriched_scalar, U_enriched_scalar);
      }
 
      // Our solution for the local enriched scalar dofs is complete. Insert into the global vector
      system.solution->insert(U_enriched_scalar, enriched_scalar_dof_indices);
    }
  }
 
  if (!global_solve)
  {
    system.solution->close();
    // Scatter solution into the current_solution which is used in error computation
    system.update();
  }
}

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::System::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::System::current_local_solution, dim, libMesh::FIFTH, DivGradExactSolution::forcing(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::index_range(), libMesh::libmesh_assert(), libMesh::DenseMatrix< T >::lu_solve(), libMesh::make_range(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::LinearImplicitSystem::reinit(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), DivGradExactSolution::scalar(), libMesh::DenseMatrix< T >::svd_solve(), libMesh::System::variable_number(), and libMesh::System::variable_type().

Referenced by assemble_divgrad(), and main().

main() [1/2]

int main

(

)

Definition at line 735 of file vector_fe_ex7.C.

{
  return 0;
}

main() [2/2]

int main	(	int	argc,
		char **	argv
	)

Definition at line 160 of file vector_fe_ex7.C.

{
  // Initialize libMesh.
  LibMeshInit init(argc, argv);
 
  // This example requires a linear solver package.
  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
                           "--enable-petsc, --enable-trilinos, or --enable-eigen");
 
  // Parse the input file.
  GetPot infile("vector_fe_ex7.in");
 
  // But allow the command line to override it.
  infile.parse_command_line(argc, argv);
 
  // Read in parameters from the command line and the input file.
  const unsigned int dimension = infile("dim", 2);
  const unsigned int grid_size = infile("grid_size", 15);
 
  // Skip higher-dimensional examples on a lower-dimensional libMesh build.
  libmesh_example_requires(dimension <= LIBMESH_DIM, dimension << "D support");
 
  // Create a mesh, with dimension to be overridden later, distributed
  // across the default MPI communicator.
  Mesh mesh(init.comm());
 
  // Use the MeshTools::Generation mesh generator to create a uniform
  // grid on the cube [-1,1]^D. To accomodate Raviart-Thomas elements, we must
  // use TRI6/7 or QUAD8/9 elements in 2d, or TET14 or HEX27 in 3d.
  const std::string elem_str = infile("element_type", std::string("TRI6"));
 
  libmesh_error_msg_if((dimension == 2 && elem_str != "TRI6" && elem_str != "TRI7" &&
                        elem_str != "QUAD8" && elem_str != "QUAD9") ||
                           (dimension == 3 && elem_str != "TET14" && elem_str != "HEX27"),
                       "You selected "
                           << elem_str
                           << " but this example must be run with TRI6, TRI7, QUAD8, or QUAD9 in 2d"
                           << " or with TET14, or HEX27 in 3d.");
 
  if (dimension == 2)
    MeshTools::Generation::build_square(
        mesh, grid_size, grid_size, -1., 1., -1., 1., Utility::string_to_enum<ElemType>(elem_str));
  else if (dimension == 3)
    MeshTools::Generation::build_cube(mesh,
                                      grid_size,
                                      grid_size,
                                      grid_size,
                                      -1.,
                                      1.,
                                      -1.,
                                      1.,
                                      -1.,
                                      1.,
                                      Utility::string_to_enum<ElemType>(elem_str));
 
  // Make sure the code is robust against nodal reorderings.
  MeshTools::Modification::permute_elements(mesh);
 
  // Create an equation systems object.
  EquationSystems equation_systems(mesh);
  const bool simplicial = (elem_str == "TRI6") || (elem_str == "TRI7") || (elem_str == "TET14");
  equation_systems.parameters.set<bool>("simplicial") = simplicial;
 
 
  // Declare the system  "DivGrad" and its variables.
  auto & system = equation_systems.add_system<System>("DivGrad");
 
  // Add the LM system
  auto & lm_system = equation_systems.add_system<LinearImplicitSystem>("Lambda");
 
  // Adds the variable "u" and "p" to "DivGrad". "u" will be our vector field
  // whereas "p" will be the scalar field.
  system.add_variable("u", FIRST, L2_RAVIART_THOMAS);
  system.add_variable("p", CONSTANT, MONOMIAL);
  // We also add a higher order version of our 'p' variable whose solution we
  // will compute using the Lagrange multiplier field and, for non-simplexes,
  // the low order 'p' solution
  system.add_variable("p_enriched", FIRST, MONOMIAL);
 
  // Add our Lagrange multiplier to the implicit system
  lm_system.add_variable("lambda", CONSTANT, SIDE_HIERARCHIC);
 
  // Give the system a pointer to the matrix assembly
  // function. This will be called when needed by the library.
  lm_system.attach_assemble_function(assemble_divgrad);
 
  // Initialize the data structures for the equation system.
  equation_systems.init();
 
  // Solve the implicit system for the Lagrange multiplier
  lm_system.solve();
 
  // Armed with our Lagrange multiplier solution, we can now compute the vector and scalar solutions
  fe_assembly(equation_systems, /*global_solve=*/false);
 
  //
  // Now we will compute our solution approximation errors
  //
 
  ExactSolution exact_sol(equation_systems);
 
  if (dimension == 2)
  {
    SolutionFunction<2> soln_func;
    SolutionGradient<2> soln_grad;
 
    // Build FunctionBase* containers to attach to the ExactSolution object.
    std::vector<FunctionBase<Number> *> sols(1, &soln_func);
    std::vector<FunctionBase<Gradient> *> grads(1, &soln_grad);
 
    exact_sol.attach_exact_values(sols);
    exact_sol.attach_exact_derivs(grads);
  }
  else if (dimension == 3)
  {
    SolutionFunction<3> soln_func;
    SolutionGradient<3> soln_grad;
 
    // Build FunctionBase* containers to attach to the ExactSolution object.
    std::vector<FunctionBase<Number> *> sols(1, &soln_func);
    std::vector<FunctionBase<Gradient> *> grads(1, &soln_grad);
 
    exact_sol.attach_exact_values(sols);
    exact_sol.attach_exact_derivs(grads);
  }
 
  // Use higher quadrature order for more accurate error results.
  int extra_error_quadrature = infile("extra_error_quadrature", 2);
  if (extra_error_quadrature)
    exact_sol.extra_quadrature_order(extra_error_quadrature);
 
  // Compute the error.
  exact_sol.compute_error("DivGrad", "u");
  exact_sol.compute_error("DivGrad", "p");
#if !defined(LIBMESH_HAVE_PETSC) || !defined(LIBMESH_USE_REAL_NUMBERS)
  if (simplicial)
#endif
    exact_sol.compute_error("DivGrad", "p_enriched");
 
  // Print out the error values.
  libMesh::out << "L2 error is: " << exact_sol.l2_error("DivGrad", "u") << std::endl;
  libMesh::out << "HDiv semi-norm error is: " << exact_sol.error_norm("DivGrad", "u", HDIV_SEMINORM)
               << std::endl;
  libMesh::out << "HDiv error is: " << exact_sol.hdiv_error("DivGrad", "u") << std::endl;
  libMesh::out << "L2 error for p is: " << exact_sol.l2_error("DivGrad", "p") << std::endl;
#if !defined(LIBMESH_HAVE_PETSC) || !defined(LIBMESH_USE_REAL_NUMBERS)
  if (simplicial)
#endif
    libMesh::out << "L2 error p_enriched is: " << exact_sol.l2_error("DivGrad", "p_enriched")
                 << std::endl;
 
#ifdef LIBMESH_HAVE_EXODUS_API
 
  // We write the file in the ExodusII format.
  ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
 
#endif // #ifdef LIBMESH_HAVE_EXODUS_API
 
  // All done.
  return 0;
}

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_divgrad(), libMesh::System::attach_assemble_function(), libMesh::ExactSolution::attach_exact_derivs(), libMesh::ExactSolution::attach_exact_values(), libMesh::MeshTools::Generation::build_cube(), libMesh::MeshTools::Generation::build_square(), libMesh::ExactSolution::compute_error(), libMesh::CONSTANT, libMesh::default_solver_package(), libMesh::ExactSolution::error_norm(), libMesh::ExactSolution::extra_quadrature_order(), fe_assembly(), libMesh::FIRST, libMesh::ExactSolution::hdiv_error(), libMesh::HDIV_SEMINORM, libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::ExactSolution::l2_error(), libMesh::L2_RAVIART_THOMAS, mesh, libMesh::MONOMIAL, libMesh::out, libMesh::EquationSystems::parameters, libMesh::MeshTools::Modification::permute_elements(), libMesh::Parameters::set(), libMesh::SIDE_HIERARCHIC, libMesh::System::solve(), and libMesh::ExodusII_IO::write_equation_systems().