Classes
class	libMesh::NumericVector< T >
	Provides a uniform interface to vector storage schemes for different linear algebra libraries. More...

class	libMesh::SparseMatrix< T >
	Generic sparse matrix. More...

Namespaces
	libMesh
	The `libMesh` namespace provides an interface to certain functionality in the library.

Functions
void	compute_residual (const NumericVector< Number > &X, NumericVector< Number > &R, NonlinearImplicitSystem &S)

void	compute_jacobian (const NumericVector< Number > &X, SparseMatrix< Number > &J, NonlinearImplicitSystem &S)

int	main (int argc, char **argv)

Function Documentation

◆ compute_jacobian()

void compute_jacobian	(	const NumericVector< Number > &	X,
		SparseMatrix< Number > &	J,
		NonlinearImplicitSystem &	S
	)

Definition at line 311 of file assembly.C.

 {
   // It is a good idea to make sure we are assembling
   // the proper system.
   libmesh_assert_equal_to(system.name(), "Poisson");
  
   // Get the DG parameters
   auto & es = system.get_equation_systems();
   const auto sigma = es.parameters.get<Real>("sigma");
   const auto epsilon = es.parameters.get<Real>("epsilon");
  
   // Get a constant reference to the mesh object.
   const MeshBase & mesh = system.get_mesh();
  
   // The dimension that we are running
   const auto dim = mesh.mesh_dimension();
  
   // 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.  We will talk more about the  DofMap
   // in future examples.
   const DofMap & dof_map = system.get_dof_map();
  
   // Get a constant reference to the Finite Element type
   // for the first (and only) variable in the system.
   FEType fe_type = dof_map.variable_type(0);
  
   // Build a Finite Element object of the specified type.
   // Note that FEVectorBase is a typedef for the templated FE
   // class.
   std::unique_ptr<FEVectorBase> fe(FEVectorBase::build(dim, fe_type));
  
   // An automatically determined Gauss quadrature rule for numerical integration.
   QGauss qrule(dim, fe_type.default_quadrature_order());
  
   // Tell the finite element object to use our quadrature rule.
   fe->attach_quadrature_rule(&qrule);
  
   // Declare a special finite element object for face integration.
   std::unique_ptr<FEVectorBase> fe_face(FEVectorBase::build(dim, fe_type));
   // And for neighbor integration
   std::unique_ptr<FEVectorBase> fe_neighbor_face(FEVectorBase::build(dim, fe_type));
  
   // Boundary integration requires one quadrature rule,
   // with dimensionality one less than the dimensionality
   // of the element.
   QGauss qface(dim - 1, fe_type.default_quadrature_order());
  
   // Tell the face finite element objects to use our
   // quadrature rule.
   fe_face->attach_quadrature_rule(&qface);
   fe_neighbor_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 auto & JxW = fe->get_JxW();
  
   // The element shape function gradients evaluated at the quadrature points.
   const auto & dphi = fe->get_dphi();
  
   // face integration points
   const auto & face_xyz = fe_face->get_xyz();
  
   // Face shape function values
   const auto & phi_face = fe_face->get_phi();
  
   // Face shape function gradients
   const auto & dphi_face = fe_face->get_dphi();
  
   // Neighbor shape function values
   const auto & phi_neighbor = fe_neighbor_face->get_phi();
  
   // Neighbor face shape function gradients
   const auto & dphi_neighbor = fe_neighbor_face->get_dphi();
  
   // face normals
   const auto & normals = fe_face->get_normals();
  
   // face JxW
   const auto & JxW_face = fe_face->get_JxW();
  
   DenseMatrix<Number> Kee;
   DenseMatrix<Number> Ken;
   DenseMatrix<Number> Kne;
   DenseMatrix<Number> Knn;
  
   // This vector 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> dof_indices_neighbor;
  
   // Now we will loop over all the elements in the mesh.
   // We will compute the element matrix and right-hand-side
   // contribution.
   //
   // Element iterators are a nice way to iterate through all the
   // elements, or all the elements that have some property.  The
   // iterator el will iterate from the first to the last element on
   // the local processor.  The iterator end_el tells us when to stop.
   // 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_elem_iterator.
   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);
  
     // Compute the element-specific data for the current
     // element.  This involves computing the location of the
     // quadrature points (q_point) and the shape functions
     // (phi, dphi) for the current element.
     fe->reinit(elem);
  
     libmesh_assert_msg(dphi.size() == dof_indices.size(),
                        "dphi size doesn't match dof_indices size");
  
     // DenseMatrix::resize() member will automatically zero out the matrix.
     Kee.resize(dof_indices.size(), dof_indices.size());
  
     // diffusion elemental jacobian
     for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
       for (std::size_t i = 0; i < dof_indices.size(); i++)
         for (std::size_t j = 0; j < dof_indices.size(); j++)
           Kee(i, j) += JxW[qp] * dphi[i][qp].contract(dphi[j][qp]);
  
     // Now we consider jacobian contributions from the sides
     for (auto side : elem->side_index_range())
     {
       // We need to compute h for penalty terms
       std::unique_ptr<const Elem> elem_side(elem->build_side_ptr(side));
       fe_face->reinit(elem, side);
       const auto elem_b_order = static_cast<unsigned int>(fe_face->get_order());
       const auto h_elem = elem->volume() / elem_side->volume() / std::pow(elem_b_order, 2);
  
       // No neighbor means we must be on a boundary
       if (!elem->neighbor_ptr(side))
       {
         for (unsigned int qp = 0; qp < qface.n_points(); qp++)
           for (std::size_t i = 0; i < dof_indices.size(); i++)
             for (std::size_t j = 0; j < dof_indices.size(); j++)
  
             {
               Kee(i, j) -= dphi_face[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
               Kee(i, j) +=
                   epsilon * phi_face[j][qp] * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
               Kee(i, j) += sigma / h_elem * phi_face[j][qp] * phi_face[i][qp] * JxW_face[qp];
             }
       }
       else // We must be on an interior side
       {
         const Elem * neighbor = elem->neighbor_ptr(side);
  
         const auto elem_id = elem->id();
         const auto neighbor_id = neighbor->id();
  
         // We don't want to erroneously add multiple contributions from the same interior face
         if ((neighbor->active() && (neighbor->level() == elem->level()) &&
              (elem_id < neighbor_id)) ||
             (neighbor->level() < elem->level()))
         {
           dof_map.dof_indices(neighbor, dof_indices_neighbor);
  
           // Make sure we have the matching quadrature points on face and neighbor
           std::vector<Point> neighbor_xyz;
           FEMap::inverse_map(elem->dim(), neighbor, face_xyz,
                              neighbor_xyz);
  
           fe_neighbor_face->reinit(neighbor, &neighbor_xyz);
  
           libmesh_assert_msg(dphi_neighbor.size() == dof_indices_neighbor.size(),
                              "dphi_neighbor size doesn't match dof_indices_neighbor size");
  
           Ken.resize(dof_indices.size(), dof_indices_neighbor.size());
           Kne.resize(dof_indices_neighbor.size(), dof_indices.size());
           Knn.resize(dof_indices_neighbor.size(), dof_indices_neighbor.size());
  
           // Now add the DG contribution to the local jacobian
           for (unsigned int qp = 0; qp < qface.n_points(); qp++)
           {
             // element-element contribution
             for (std::size_t i = 0; i < dof_indices.size(); i++)
               for (std::size_t j = 0; j < dof_indices.size(); j++)
               {
                 Kee(i, j) -= 0.5 * dphi_face[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
                 Kee(i, j) +=
                     epsilon * 0.5 * phi_face[j][qp] * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
                 Kee(i, j) += sigma / h_elem * phi_face[j][qp] * phi_face[i][qp] * JxW_face[qp];
               }
             // element-neighbor contribution
             for (std::size_t i = 0; i < dof_indices.size(); i++)
               for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
               {
                 Ken(i, j) -=
                     0.5 * dphi_neighbor[j][qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
                 Ken(i, j) += epsilon * 0.5 * -phi_neighbor[j][qp] * dphi_face[i][qp] * normals[qp] *
                              JxW_face[qp];
                 Ken(i, j) += sigma / h_elem * -phi_neighbor[j][qp] * phi_face[i][qp] * JxW_face[qp];
               }
             // Neighbor-element contribution
             for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
               for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
               {
                 Kne(i, j) +=
                     0.5 * dphi_face[j][qp] * normals[qp] * phi_neighbor[i][qp] * JxW_face[qp];
                 Kne(i, j) -= epsilon * 0.5 * phi_face[j][qp] * dphi_neighbor[i][qp] * normals[qp] *
                              JxW_face[qp];
                 Kne(i, j) -= sigma / h_elem * phi_face[j][qp] * phi_neighbor[i][qp] * JxW_face[qp];
               }
             // Neighbor-neighbor contribution
             for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
               for (std::size_t j = 0; j < dof_indices_neighbor.size(); j++)
               {
                 Knn(i, j) +=
                     0.5 * dphi_neighbor[j][qp] * normals[qp] * phi_neighbor[i][qp] * JxW_face[qp];
                 Knn(i, j) -= epsilon * 0.5 * -phi_neighbor[j][qp] * dphi_neighbor[i][qp] *
                              normals[qp] * JxW_face[qp];
                 Knn(i, j) -=
                     sigma / h_elem * -phi_neighbor[j][qp] * phi_neighbor[i][qp] * JxW_face[qp];
               }
           }
  
           J.add_matrix(Ken, dof_indices, dof_indices_neighbor);
           J.add_matrix(Kne, dof_indices_neighbor, dof_indices);
           J.add_matrix(Knn, dof_indices_neighbor, dof_indices_neighbor);
         } // whether we've done this internal side integral before
       }   // whether we're on an internal side
     }     // loop over sides
  
     J.add_matrix(Kee, dof_indices, dof_indices);
   } // element loop
 }

References libMesh::Elem::active(), libMesh::MeshBase::active_local_element_ptr_range(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::FEMap::inverse_map(), libMesh::Elem::level(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::System::name(), libMesh::Elem::neighbor_ptr(), libMesh::EquationSystems::parameters, std::pow(), libMesh::Real, libMesh::DenseMatrix< T >::resize(), and libMesh::DofMap::variable_type().

Referenced by libMesh::FirstOrderUnsteadySolver::compute_second_order_eqns(), PoissonSystem::element_time_derivative(), LaplaceSystem::element_time_derivative(), HeatSystem::element_time_derivative(), main(), and LaplaceSystem::side_constraint().

◆ compute_residual()

void compute_residual	(	const NumericVector< Number > &	X,
		NumericVector< Number > &	R,
		NonlinearImplicitSystem &	S
	)

Vectors to hold the local solution degree of freedom values

Vector to hold the local solution

Definition at line 26 of file assembly.C.

 {
   // It is a good idea to make sure we are assembling
   // the proper system.
   libmesh_assert_equal_to(system.name(), "Poisson");
  
   // Get the DG parameters
   auto & es = system.get_equation_systems();
   const auto sigma = es.parameters.get<Real>("sigma");
   const auto epsilon = es.parameters.get<Real>("epsilon");
  
   // Get a constant reference to the mesh object.
   const MeshBase & mesh = system.get_mesh();
  
   // The dimension that we are running
   const auto dim = mesh.mesh_dimension();
  
   // 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.  We will talk more about the  DofMap
   // in future examples.
   const DofMap & dof_map = system.get_dof_map();
  
   // Get a constant reference to the Finite Element type
   // for the first (and only) variable in the system.
   FEType fe_type = dof_map.variable_type(0);
  
   // Build a Finite Element object of the specified type.
   // Note that FEVectorBase is a typedef for the templated FE
   // class.
   std::unique_ptr<FEVectorBase> fe(FEVectorBase::build(dim, fe_type));
  
   // An automatically determined Gauss quadrature rule for numerical integration.
   QGauss qrule(dim, fe_type.default_quadrature_order());
  
   // Tell the finite element object to use our quadrature rule.
   fe->attach_quadrature_rule(&qrule);
  
   // Declare a special finite element object for face integration.
   std::unique_ptr<FEVectorBase> fe_face(FEVectorBase::build(dim, fe_type));
   // And for neighbor integration
   std::unique_ptr<FEVectorBase> fe_neighbor_face(FEVectorBase::build(dim, fe_type));
  
   // Boundary integration requires one quadrature rule,
   // with dimensionality one less than the dimensionality
   // of the element.
   QGauss qface(dim - 1, fe_type.default_quadrature_order());
  
   // Tell the face finite element objects to use our
   // quadrature rule.
   fe_face->attach_quadrature_rule(&qface);
   fe_neighbor_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 auto & JxW = fe->get_JxW();
  
   // element integration points
   const auto & xyz = fe->get_xyz();
  
   // The element shape function values evaluated at the quadrature points
   const auto & phi = fe->get_phi();
  
   // The element shape function gradients evaluated at the quadrature points.
   const auto & dphi = fe->get_dphi();
  
   // face integration points
   const auto & face_xyz = fe_face->get_xyz();
  
   // Face shape function values
   const auto & phi_face = fe_face->get_phi();
  
   // Face shape function gradients
   const auto & dphi_face = fe_face->get_dphi();
  
   // Neighbor shape function values
   const auto & phi_neighbor = fe_neighbor_face->get_phi();
  
   // Neighbor face shape function gradients
   const auto & dphi_neighbor = fe_neighbor_face->get_dphi();
  
   // face normals
   const auto & normals = fe_face->get_normals();
  
   // face JxW
   const auto & JxW_face = fe_face->get_JxW();
  
   DenseVector<Number> Fe;
   DenseVector<Number> Fn;
  
   // This vector 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> dof_indices_neighbor;
  
   std::vector<Number> dof_u;
   std::vector<Number> dof_u_neighbor;
  
   std::vector<VectorValue<Number>> u;
   std::vector<TensorValue<Number>> grad_u;
   std::vector<VectorValue<Number>> u_neighbor;
   std::vector<TensorValue<Number>> grad_u_neighbor;
  
   // Now we will loop over all the elements in the mesh.
   // We will compute the element matrix and right-hand-side
   // contribution.
   //
   // Element iterators are a nice way to iterate through all the
   // elements, or all the elements that have some property.  The
   // iterator el will iterate from the first to the last element on
   // the local processor.  The iterator end_el tells us when to stop.
   // 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_elem_iterator.
   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);
  
     // Compute the element-specific data for the current
     // element.  This involves computing the location of the
     // quadrature points (q_point) and the shape functions
     // (phi, dphi) for the current element.
     fe->reinit(elem);
  
     libmesh_assert_msg(dphi.size() == dof_indices.size(),
                        "dphi size doesn't match dof_indices size");
  
     // DenseVector::resize() member will automatically zero out the vector.
     Fe.resize(dof_indices.size());
  
     // Get the local solution vector
     X.get(dof_indices, dof_u);
  
     // build the element solution gradient
     grad_u.resize(qrule.n_points());
     for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
     {
       grad_u[qp] = 0;
       for (std::size_t i = 0; i < dof_indices.size(); i++)
         grad_u[qp] += dof_u[i] * dphi[i][qp];
     }
  
     for (unsigned int qp = 0; qp < qrule.n_points(); qp++)
     {
       auto forcing_func = VectorValue<Number>(forcing_function(0, xyz[qp](0), xyz[qp](1)),
                                               forcing_function(1, xyz[qp](0), xyz[qp](1)),
                                               0);
       for (std::size_t i = 0; i < dof_indices.size(); i++)
       {
         // diffusion elemental residual
         Fe(i) += JxW[qp] * dphi[i][qp].contract(grad_u[qp]);
  
         // forcing function
         Fe(i) += JxW[qp] * phi[i][qp] * forcing_func;
       }
     }
  
     // Now we consider residual contributions from the sides
     for (auto side : elem->side_index_range())
     {
       // We need to compute h for penalty terms
       std::unique_ptr<const Elem> elem_side(elem->build_side_ptr(side));
       fe_face->reinit(elem, side);
       const auto elem_b_order = static_cast<unsigned int>(fe_face->get_order());
       const auto h_elem = elem->volume() / elem_side->volume() / std::pow(elem_b_order, 2);
  
       // build the face solution value and gradient
       u.resize(qface.n_points());
       grad_u.resize(qface.n_points());
       for (unsigned int qp = 0; qp < qface.n_points(); qp++)
       {
         u[qp] = 0;
         grad_u[qp] = 0;
         for (std::size_t i = 0; i < dof_indices.size(); i++)
         {
           u[qp] += dof_u[i] * phi_face[i][qp];
           grad_u[qp] += dof_u[i] * dphi_face[i][qp];
         }
       }
  
       // No neighbor means we must be on a boundary
       if (!elem->neighbor_ptr(side))
       {
         for (unsigned int qp = 0; qp < qface.n_points(); qp++)
         {
           auto fn = VectorValue<Number>(exact_solution(0, face_xyz[qp](0), face_xyz[qp](1)),
                                         exact_solution(1, face_xyz[qp](0), face_xyz[qp](1)),
                                         0);
           for (std::size_t i = 0; i < dof_indices.size(); i++)
           {
             Fe(i) -= grad_u[qp] * normals[qp] * phi_face[i][qp] * JxW_face[qp];
             Fe(i) += epsilon * (u[qp] - fn) * dphi_face[i][qp] * normals[qp] * JxW_face[qp];
             Fe(i) += sigma / h_elem * (u[qp] - fn) * phi_face[i][qp] * JxW_face[qp];
           }
         }
       }
       else // We must be on an interior side
       {
         const Elem * neighbor = elem->neighbor_ptr(side);
  
         const auto elem_id = elem->id();
         const auto neighbor_id = neighbor->id();
  
         // We don't want to erroneously add multiple contributions from the same interior face
         if ((neighbor->active() && (neighbor->level() == elem->level()) &&
              (elem_id < neighbor_id)) ||
             (neighbor->level() < elem->level()))
         {
           dof_map.dof_indices(neighbor, dof_indices_neighbor);
  
           // Make sure we have the matching quadrature points on face and neighbor
           std::vector<Point> neighbor_xyz;
           FEMap::inverse_map(elem->dim(), neighbor, face_xyz,
                              neighbor_xyz);
  
           fe_neighbor_face->reinit(neighbor, &neighbor_xyz);
  
           libmesh_assert_msg(dphi_neighbor.size() == dof_indices_neighbor.size(),
                              "dphi_neighbor size doesn't match dof_indices_neighbor size");
  
           Fn.resize(dof_indices_neighbor.size());
           X.get(dof_indices_neighbor, dof_u_neighbor);
  
           // build the neighbor solution value and gradient
           u_neighbor.resize(qface.n_points());
           grad_u_neighbor.resize(qface.n_points());
           for (unsigned int qp = 0; qp < qface.n_points(); qp++)
           {
             u_neighbor[qp] = 0;
             grad_u_neighbor[qp] = 0;
             for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
             {
               u_neighbor[qp] += dof_u_neighbor[i] * phi_neighbor[i][qp];
               grad_u_neighbor[qp] += dof_u_neighbor[i] * dphi_neighbor[i][qp];
             }
           }
  
           // Now add the DG contribution to the local residual
           for (unsigned int qp = 0; qp < qface.n_points(); qp++)
           {
             // element contribution
             for (std::size_t i = 0; i < dof_indices.size(); i++)
             {
               Fe(i) -= 0.5 * (grad_u[qp] * normals[qp] + grad_u_neighbor[qp] * normals[qp]) *
                        phi_face[i][qp] * JxW_face[qp];
               Fe(i) += epsilon * 0.5 * (u[qp] - u_neighbor[qp]) * dphi_face[i][qp] * normals[qp] *
                        JxW_face[qp];
               Fe(i) += sigma / h_elem * (u[qp] - u_neighbor[qp]) * phi_face[i][qp] * JxW_face[qp];
             }
             // Neighbor contribution
             for (std::size_t i = 0; i < dof_indices_neighbor.size(); i++)
             {
               Fn(i) += 0.5 * (grad_u[qp] * normals[qp] + grad_u_neighbor[qp] * normals[qp]) *
                        phi_neighbor[i][qp] * JxW_face[qp];
               Fn(i) -= epsilon * 0.5 * (u[qp] - u_neighbor[qp]) * dphi_neighbor[i][qp] *
                        normals[qp] * JxW_face[qp];
               Fn(i) -=
                   sigma / h_elem * (u[qp] - u_neighbor[qp]) * phi_neighbor[i][qp] * JxW_face[qp];
             }
           }
  
           R.add_vector(Fn, dof_indices_neighbor);
         } // whether we've done this internal side integral before
       }   // whether we're on an internal side
     }     // loop over sides
  
     R.add_vector(Fe, dof_indices);
   } // element loop
 }

References libMesh::Elem::active(), libMesh::MeshBase::active_local_element_ptr_range(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), exact_solution(), forcing_function(), libMesh::Parameters::get(), libMesh::NumericVector< T >::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::FEMap::inverse_map(), libMesh::Elem::level(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::System::name(), libMesh::Elem::neighbor_ptr(), libMesh::EquationSystems::parameters, std::pow(), libMesh::Real, libMesh::DenseVector< T >::resize(), and libMesh::DofMap::variable_type().

Referenced by main().

◆ main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 62 of file vector_fe_ex5.C.

 {
   // Initialize libraries.
   LibMeshInit init(argc, argv);
  
   // This example requires a non-linear solver package.
   libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS ||
                                libMesh::default_solver_package() == TRILINOS_SOLVERS,
                            "--enable-petsc or --enable-trilinos");
  
   // Parse input file
   GetPot input_file("vector_fe_ex5.in");
  
   // Read DG parameters
   const Real epsilon = input_file("epsilon", -1);
   const Real sigma = input_file("sigma", 6);
  
   // Read mesh size
   const std::size_t nx = input_file("nx", 15);
   const std::size_t ny = input_file("ny", 15);
  
   // Brief message to the user regarding the program name
   // and command line arguments.
   libMesh::out << "Running " << argv[0];
  
   for (int i = 1; i < argc; i++)
     libMesh::out << " " << argv[i];
  
   libMesh::out << std::endl << std::endl;
  
   // Skip this 2D example if libMesh was compiled as 1D-only.
   libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
  
   // Create a mesh, with dimension to be overridden later, on the
   // default MPI communicator.
   Mesh mesh(init.comm());
  
   // Use the MeshTools::Generation mesh generator to create a uniform
   // 2D grid on the square [-1,1]^2.  We instruct the mesh generator
   // to build a mesh of 15x15 QUAD4 elements. Note that at the end of this
   // fuction call, the mesh will be prepared and remote elements will be deleted
   MeshTools::Generation::build_square(mesh, nx, ny, -1., 1., -1., 1., QUAD4);
  
   // Print information about the mesh to the screen.
   mesh.print_info();
  
   // Create an equation systems object.
   EquationSystems equation_systems(mesh);
   equation_systems.parameters.set<Real>("epsilon") = epsilon;
   equation_systems.parameters.set<Real>("sigma") = sigma;
  
   // Declare the Poisson system and its variables.
   // The Poisson system is another example of a steady system.
   auto & nl_system = equation_systems.add_system<NonlinearImplicitSystem>("Poisson");
  
   // Adds the variable "u" to "Poisson".  "u"
   // will be approximated using first-order vector monomial elements.
   // Since the mesh is 2-D, "u" will have two components.
   nl_system.add_variable("u", FIRST, MONOMIAL_VEC);
  
   // Set the residual and Jacobian evaluation functions
   auto & nl_solver = *nl_system.nonlinear_solver;
   nl_solver.residual = compute_residual;
   nl_solver.jacobian = compute_jacobian;
  
   // Initialize the data structures for the equation system.
   equation_systems.init();
  
   // Prints information about the system to the screen.
   equation_systems.print_info();
  
   nl_system.solve();
  
 #ifdef LIBMESH_HAVE_EXODUS_API
   ExodusII_IO(mesh).write_equation_systems("out.e", equation_systems);
 #endif
  
   // All done.
   return 0;
 }

References libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), libMesh::MeshTools::Generation::build_square(), compute_jacobian(), compute_residual(), libMesh::default_solver_package(), libMesh::FIRST, libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), mesh, libMesh::MONOMIAL_VEC, libMesh::out, libMesh::EquationSystems::parameters, libMesh::PETSC_SOLVERS, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD4, libMesh::Real, libMesh::Parameters::set(), libMesh::TRILINOS_SOLVERS, and libMesh::MeshOutput< MT >::write_equation_systems().

Classes

Namespaces

Functions

Function Documentation

◆ compute_jacobian()

◆ compute_residual()

◆ main()