vector_fe_ex6.C File Reference

Go to the source code of this file.

Functions
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))

Function Documentation

assemble_divgrad() [1/2]

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

Referenced by main().

assemble_divgrad() [2/2]

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

Definition at line 260 of file vector_fe_ex6.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, DivGradExactSolution::forcing(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::ImplicitSystem::get_system_matrix(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, DivGradExactSolution::scalar(), and libMesh::System::variable_number().

{

  // It is a good idea to make sure we are assembling
  // the proper system.
  libmesh_assert_equal_to (system_name, "DivGrad");

  // Retrieve our boundary condition type. If not Neumann, then it is Dirichlet.
  const bool neumann = es.parameters.get<bool>("neumann");

  // Get a constant reference to the mesh object.
  const MeshBase & mesh = es.get_mesh();

  // The dimension that we are running.
  const unsigned int dim = mesh.mesh_dimension();

  // Get a reference to the LinearImplicitSystem we are solving.
  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("DivGrad");

  // A reference to the  DofMap object for this system.  The  DofMap
  // object handles the index translation from node and element numbers
  // to degree of freedom numbers.
  const DofMap & dof_map = system.get_dof_map();

  // Get a constant reference to the Finite Element type
  // for the two variables in the system.
  FEType vector_fe_type = dof_map.variable_type(system.variable_number("u"));
  FEType scalar_fe_type = dof_map.variable_type(system.variable_number("p"));

  // Build two Finite Element objects, one of each specified type. Since the
  // FEBase::build() member dynamically creates memory we will
  // store the object as a std::unique_ptr<FEBase>. This can be thought
  // of as a pointer that will clean up after itself. Introduction Example 4
  // describes some advantages of  std::unique_ptr's in the context of
  // quadrature rules.
  std::unique_ptr<FEVectorBase> vector_fe (FEVectorBase::build(dim, vector_fe_type));
  std::unique_ptr<FEBase> scalar_fe (FEBase::build(dim, scalar_fe_type));

  // A just-high-enough Gauss quadrature rule for numerical integration.
  QGauss qrule (dim, vector_fe_type.default_quadrature_order());

  // Tell the finite element objects to use our quadrature rule.
  vector_fe->attach_quadrature_rule (&qrule);
  scalar_fe->attach_quadrature_rule (&qrule);

  // Declare a special finite element object for boundary integration.
  std::unique_ptr<FEVectorBase> vector_fe_face (FEVectorBase::build(dim, vector_fe_type));

  // Boundary integration requires one quadrature rule with dimensionality one
  // less than the dimensionality of the element.
  QGauss qface(dim-1, vector_fe_type.default_quadrature_order());

  // Tell the finite element object to use our quadrature rule.
  vector_fe_face->attach_quadrature_rule (&qface);

  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.
  //
  // The element Jacobian * quadrature weight at each integration point.
  const std::vector<Real> & JxW = vector_fe->get_JxW();

  // The physical XY locations of the quadrature points on the element.
  // These might be useful for evaluating spatially varying material
  // properties at the quadrature points.
  const std::vector<Point> & q_point = vector_fe->get_xyz();

  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<RealGradient>> & vector_phi = vector_fe->get_phi();
  const std::vector<std::vector<Real>> & scalar_phi = scalar_fe->get_phi();

  // The divergence of the element vector shape functions evaluated at the
  // quadrature points.
  const std::vector<std::vector<Real>> & div_vector_phi = vector_fe->get_div_phi();

  // Define data structures to contain the element matrix
  // and right-hand-side vector contribution.  Following
  // basic finite element terminology we will denote these
  // "Ke" and "Fe".  These datatypes are templated on
  //  Number, which allows the same code to work for real
  // or complex numbers.
  DenseMatrix<Number> Ke;
  DenseVector<Number> Fe;

  // These vectors will hold the degree of freedom indices for
  // the element.  These define where in the global system
  // the element degrees of freedom get mapped.
  std::vector<dof_id_type> dof_indices;
  std::vector<dof_id_type> vector_dof_indices;
  std::vector<dof_id_type> scalar_dof_indices;
  std::vector<dof_id_type> lambda_dof_indices;

  // The global system matrix
  SparseMatrix<Number> & matrix = system.get_system_matrix();

  // Now we will loop over all the elements in the mesh.
  // We will compute the element matrix and right-hand-side
  // contribution.
  //
  // Element ranges are a nice way to iterate through all the
  // elements, or all the elements that have some property.  The
  // range will iterate from the first to the last element on
  // the local processor.
  // It is smart to make this one const so that we don't accidentally
  // mess it up!  In case users later modify this program to include
  // refinement, we will be safe and will only consider the active
  // elements; hence we use a variant of the
  // active_local_element_ptr_range.
  for (const auto & elem : mesh.active_local_element_ptr_range())
    {
      // Get the degree of freedom indices for the
      // current element.  These define where in the global
      // matrix and right-hand-side this element will
      // contribute to.
      dof_map.dof_indices (elem, dof_indices);
      dof_map.dof_indices (elem, vector_dof_indices, system.variable_number("u"));
      dof_map.dof_indices (elem, scalar_dof_indices, system.variable_number("p"));
      if (neumann)
        dof_map.dof_indices (elem, lambda_dof_indices, system.variable_number("l"));

      // Cache the number of degrees of freedom, in total and for each
      // variable, on this element, for use as array and loop bounds later.
      // We use cast_int to explicitly convert from size() (which may be
      // 64-bit) to unsigned int (which may be 32-bit but which is definitely
      // enough to count *local* degrees of freedom.
      const unsigned int n_dofs =
        cast_int<unsigned int>(dof_indices.size());
      const unsigned int vector_n_dofs =
        cast_int<unsigned int>(vector_dof_indices.size());
      const unsigned int scalar_n_dofs =
        cast_int<unsigned int>(scalar_dof_indices.size());
      const unsigned int lambda_n_dofs =
        cast_int<unsigned int>(lambda_dof_indices.size());

      // Compute the element-specific data for the current
      // element.  This involves computing the location of the
      // quadrature points (q_point) and the shape functions
      // and their divergences for the current element.
      vector_fe->reinit (elem);
      scalar_fe->reinit (elem);

      // The total number of degrees of freedom is just the sum of the number
      // of degrees of freedom per variable. We should also have the same
      // number of degrees of freedom as shape functions for each variable.
      libmesh_assert_equal_to (n_dofs, vector_n_dofs + scalar_n_dofs + lambda_n_dofs);
      libmesh_assert_equal_to (vector_n_dofs, vector_phi.size());
      libmesh_assert_equal_to (scalar_n_dofs, scalar_phi.size());

      // Zero the element matrix and right-hand side before
      // summing them.  We use the resize member here because
      // the number of degrees of freedom might have changed from
      // the last element.  Note that this will be the case if the
      // element type is different (i.e. the last element was a
      // triangle, now we are on a quadrilateral).

      // The  DenseMatrix::resize() and the  DenseVector::resize()
      // members will automatically zero out the matrix  and vector.
      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);

      // Now loop over the quadrature points.  This handles
      // the numeric integration.
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {

          // Now we will build the element matrix.
          // The upper-left block involves a double loop to integrate the
          // vector test functions (i) against the vector trial functions (j).
          for (unsigned int i = 0; i != vector_n_dofs; i++)
            for (unsigned int j = 0; j != vector_n_dofs; j++)
              {
                Ke(i, j) += JxW[qp]*(vector_phi[i][qp]*vector_phi[j][qp]);
              }

          // The upper-right block involves a double loop to integrate the
          // divergence of the vector test functions (i) against the scalar
          // trial functions (l).
          for (unsigned int i = 0; i != vector_n_dofs; i++)
            for (unsigned int l = 0; l != scalar_n_dofs; l++)
              {
                Ke(i, l + vector_n_dofs) -= JxW[qp]*(div_vector_phi[i][qp]*scalar_phi[l][qp]);
              }

          // The lower-left block involves a double loop to integrate the
          // scalar test functions (k) against the divergence of the vector
          // trial functions (j).
          for (unsigned int k = 0; k != scalar_n_dofs; k++)
            for (unsigned int j = 0; j != vector_n_dofs; j++)
              {
                Ke(k + vector_n_dofs, j) -= JxW[qp]*(div_vector_phi[j][qp]*scalar_phi[k][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).
          {
            // The location of the current quadrature point.
            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);

            // Loop to integrate the scalar test functions (k) against the
            // forcing function.
            for (unsigned int k = 0; k != scalar_n_dofs; k++)
              {
                Fe(k + vector_n_dofs) -= JxW[qp]*f*scalar_phi[k][qp];
              }
          }

          // We have now reached the end of the RHS summation. In addition,
          // however, since the scalar variable is defined only up
          // to an additive constant with purely Neumann boundary conditions, we
          // constrain the integral of the scalar variable to the integral
          // of the exact solution we seek.
          {
            // The location of the current quadrature point.
            const Real x = q_point[qp](0);
            const Real y = q_point[qp](1);
            const Real z = q_point[qp](2);

            // The value of the scalar variable.
            Real scalar_value = 0;
            if (dim == 2)
              scalar_value = DivGradExactSolution().scalar(x, y);
            else if (dim == 3)
              scalar_value = DivGradExactSolution().scalar(x, y, z);

            // A double loop to integrate the
            // scalar test functions (k) against the Lagrange dof (n).
            for (unsigned int k = 0; k != scalar_n_dofs; k++)
              for (unsigned int n = 0; n != lambda_n_dofs; n++)
                {
                  Ke(k + vector_n_dofs, n + vector_n_dofs + scalar_n_dofs) += JxW[qp]*scalar_phi[k][qp];
                }

            // A double loop to integrate the Lagrange dof (m) against the
            // scalar trial functions (l).
            for (unsigned int m = 0; m != lambda_n_dofs; m++)
              for (unsigned int l = 0; l != scalar_n_dofs; l++)
                {
                  Ke(m + vector_n_dofs + scalar_n_dofs, l + vector_n_dofs) += JxW[qp]*scalar_phi[l][qp];
                }

            // Loop to integrate the exact solution for the scalar variable.
            for (unsigned int m = 0; m != lambda_n_dofs; m++)
              {
                Fe(m + vector_n_dofs + scalar_n_dofs) += JxW[qp]*scalar_value;
              }
          }
        }

      // We have now reached the end of the quadrature point loop, so
      // the interior element integration has
      // been completed.  However, we have not yet addressed
      // boundary conditions.  For this Poisson example, we consider either
      // Dirichlet or Neumann for the scalar solution field p. Note that in
      // the mixed formulation a Neumann condition for p corresponds to a
      // Dirichlet condition for u.
      {

        // The following loop is over the sides of the element.
        // If the element has no neighbor on a side then that
        // side MUST live on a boundary of the domain.
        for (auto side : elem->side_index_range())
          if (elem->neighbor_ptr(side) == nullptr)
            {
              // The value of the shape functions at the quadrature points.
              const std::vector<std::vector<RealGradient>> & vector_phi_face = vector_fe_face->get_phi();

              // The Jacobian * Quadrature Weight at the quadrature
              // points on the face.
              const std::vector<Real> & JxW_face = vector_fe_face->get_JxW();

              // The XYZ locations (in physical space) of, and the normals at,
              // the quadrature points on the face.  This is where
              // we will interpolate the boundary value function.
              const std::vector<Point> & qface_point = vector_fe_face->get_xyz();
              const std::vector<Point> & normals = vector_fe_face->get_normals();

              // Compute the vector shape function values on the element face.
              vector_fe_face->reinit(elem, side);

              // Some shape functions will be 0 on the face, but for ease of
              // indexing and generality of code we loop over them anyway.
              libmesh_assert_equal_to (vector_n_dofs, vector_phi_face.size());

              // Loop over the face quadrature points for integration.
              for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  // The location on the boundary of the current
                  // face quadrature point.
                  const Real xf = qface_point[qp](0);
                  const Real yf = qface_point[qp](1);
                  const Real zf = qface_point[qp](2);

                  if (neumann)
                    {
                      // The boundary value for the vector variable.
                      RealGradient vector_value;
                      if (dim == 2)
                        vector_value = DivGradExactSolution()(xf, yf);
                      else if (dim == 3)
                        vector_value = DivGradExactSolution()(xf, yf, zf);

                      // We use the penalty method to set the flux of the vector
                      // variable at the boundary, i.e. the RT vector boundary dof.
                      const Real penalty = 1.e10;

                      // A double loop to integrate the normal component of the
                      // vector test functions (i) against the normal component of
                      // the vector trial functions (j).
                      for (unsigned int i = 0; i != vector_n_dofs; i++)
                        for (unsigned int j = 0; j != vector_n_dofs; j++)
                          {
                            Ke(i, j) += JxW_face[qp]*penalty*vector_phi_face[i][qp]*
                                        normals[qp]*vector_phi_face[j][qp]*normals[qp];
                          }

                      // Loop to integrate the normal component of the vector test
                      // functions (i) against the normal component of the
                      // exact solution for the vector variable.
                      for (unsigned int i = 0; i != vector_n_dofs; i++)
                        {
                          Fe(i) += JxW_face[qp]*penalty*vector_phi_face[i][qp]*normals[qp]*
                                   vector_value*normals[qp];
                        }
                    }
                  else
                    {
                      // 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);

                      // Loop to integrate the normal component of the vector test
                      // functions (i) against the exact solution for the scalar variable.
                      for (unsigned int i = 0; i != vector_n_dofs; i++)
                        {
                          Fe(i) += -JxW_face[qp]*vector_phi_face[i][qp]*normals[qp]*scalar_value;
                        }
                    }
                }
            }
      }

      // We have now finished the quadrature point loop,
      // and have therefore applied all the boundary conditions.

      // If this assembly program were to be used on an adaptive mesh,
      // we would have to apply any hanging node constraint equations.
      dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);

      // The element matrix and right-hand-side are now built
      // for this element.  Add them to the global matrix and
      // right-hand-side vector.  The  SparseMatrix::add_matrix()
      // and  NumericVector::add_vector() members do this for us.
      matrix.add_matrix (Ke, dof_indices);
      system.rhs->add_vector (Fe, dof_indices);
    }

  // All done!
}

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 83 of file vector_fe_ex6.C.

References libMesh::EquationSystems::add_system(), assemble_divgrad(), 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(), libMesh::FIFTH, libMesh::FIRST, libMesh::ExactSolution::hdiv_error(), libMesh::HDIV_SEMINORM, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::ExactSolution::l2_error(), libMesh::L2_HIERARCHIC, mesh, libMesh::MONOMIAL, libMesh::out, libMesh::EquationSystems::parameters, libMesh::MeshTools::Modification::permute_elements(), libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::RAVIART_THOMAS, libMesh::SCALAR, libMesh::Parameters::set(), and libMesh::ExodusII_IO::write_equation_systems().

{
  // 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_ex6.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.");

  const std::string bc_str = infile("boundary_condition", std::string("neumann"));
  libmesh_error_msg_if(
      (bc_str != "neumann") && (bc_str != "dirichlet"),
      "You selected '" << bc_str << "', however, the valid options are 'dirichlet' or 'neumann'");
  const bool neumann = (bc_str == "neumann");

  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);

  // Print information about the mesh to the screen.
  mesh.print_info();

  // Create an equation systems object.
  EquationSystems equation_systems (mesh);
  equation_systems.parameters.set<bool>("neumann") = neumann;

  // Declare the system  "DivGrad" and its variables.
  LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem>("DivGrad");

  // Set the FE approximation order for the vector and scalar field variables.
  const Order vector_order = static_cast<Order>(infile("order", 1u));
  const Order scalar_order = static_cast<Order>(vector_order - 1u);

  libmesh_error_msg_if(vector_order < FIRST || vector_order > ((dimension == 3) ? FIRST : FIFTH),
                       "You selected: " << vector_order <<
                       " but this example must be run with either 1 <= order <= 5 in 2d"
                       " or with order 1 in 3d.");

  // Adds the variables "u" and "p" to "DivGrad". "u" will be our vector field
  // whereas "p" will be the scalar field.
  system.add_variable("u", vector_order, RAVIART_THOMAS);
  system.add_variable("p", scalar_order, scalar_order == CONSTANT ? MONOMIAL : L2_HIERARCHIC);

  // Add a scalar Lagrange multiplier to remove the nullspace if imposing the Neumann condition.
  if (neumann)
    system.add_variable("l", FIRST, SCALAR);

  // Give the system a pointer to the matrix assembly
  // function. This will be called when needed by the library.
  system.attach_assemble_function(assemble_divgrad);

  // Initialize the data structures for the equation system.
  equation_systems.init();

  // Prints information about the system to the screen.
  equation_systems.print_info();

  // Solve the system "DivGrad". Note that calling this
  // member will assemble the linear system and invoke
  // the default numerical solver.
  system.solve();

  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);
  exact_sol.extra_quadrature_order(extra_error_quadrature);

  // Compute the error.
  exact_sol.compute_error("DivGrad", "u");
  exact_sol.compute_error("DivGrad", "p");

  // Print out the error values.
  libMesh::out << "~~ Vector field (u) ~~"
               << std::endl;
  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 << "~~ Scalar field (p) ~~"
               << std::endl;
  libMesh::out << "L2 error is: "
               << exact_sol.l2_error("DivGrad", "p")
               << 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;
}