LaplaceYoung Class Reference

A class which provides the residual and jacobian assembly functions for the Laplace-Young system of equations. More...

Inheritance diagram for LaplaceYoung:

Public Member Functions
	LaplaceYoung ()
virtual void	residual (const NumericVector< Number > &X, NumericVector< Number > &R, NonlinearImplicitSystem &S)
	Function which computes the residual. More...
virtual void	jacobian (const NumericVector< Number > &X, SparseMatrix< Number > &J, NonlinearImplicitSystem &S)
	Function which computes the jacobian. More...
virtual void	postcheck (const NumericVector< Number > &old_soln, NumericVector< Number > &search_direction, NumericVector< Number > &new_soln, bool &changed_search_direction, bool &changed_new_soln, NonlinearImplicitSystem &S)
	Function which performs a postcheck on the solution. More...
virtual void	residual (const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S)=0
	Residual function. More...
virtual void	jacobian (const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)=0
	Jacobian function. More...
virtual void	postcheck (const NumericVector< Number > &old_soln, NumericVector< Number > &search_direction, NumericVector< Number > &new_soln, bool &changed_search_direction, bool &changed_new_soln, sys_type &S)=0
	This interface, which is inspired by PETSc's, passes the user: .) A constant reference to the "old" solution vector (previous Newton iterate), .) A pointer to the search direction (update) vector, and .) The "proposed" new solution vector. More...

Private Attributes
Real	_kappa
Real	_sigma
Real	_gamma

Detailed Description

A class which provides the residual and jacobian assembly functions for the Laplace-Young system of equations.

Definition at line 77 of file miscellaneous_ex3.C.

Constructor & Destructor Documentation

LaplaceYoung()

LaplaceYoung::LaplaceYoung ( )

inline

Definition at line 83 of file miscellaneous_ex3.C.

                 :
    _kappa(1.),
    _sigma(0.2),
    _gamma(1.0)
  {}

Member Function Documentation

jacobian() [1/2]

void LaplaceYoung::jacobian ( const NumericVector< Number > &  X, SparseMatrix< Number > &  J, NonlinearImplicitSystem &  S  )

virtual

Function which computes the jacobian.

Definition at line 464 of file miscellaneous_ex3.C.

{
  // Get a reference to the equation system.
  EquationSystems & es = sys.get_equation_systems();
 
  // 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 NonlinearImplicitSystem we are solving
  NonlinearImplicitSystem & system =
    es.get_system<NonlinearImplicitSystem>("Laplace-Young");
 
  // 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.  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.
  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
 
  // A 5th order Gauss quadrature rule for numerical integration.
  QGauss qrule (dim, FIFTH);
 
  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule (&qrule);
 
  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.
  // We begin with the element Jacobian * quadrature weight at each
  // integration point.
  const std::vector<Real> & JxW = fe->get_JxW();
 
  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<Real>> & phi = fe->get_phi();
 
  // The element shape function gradients evaluated at the quadrature
  // points.
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
 
  // Define data structures to contain the Jacobian element matrix.
  // Following basic finite element terminology we will denote these
  // "Ke". More detail is in example 3.
  DenseMatrix<Number> Ke;
 
  // 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;
 
  // Now we will loop over all the active elements in the mesh which
  // are local to this processor.
  // We will compute the element Jacobian contribution.
  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);
 
      // Zero the element Jacobian 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).
      const unsigned int n_dofs =
        cast_int<unsigned int>(dof_indices.size());
      Ke.resize (n_dofs, n_dofs);
 
      // Now we will build the element Jacobian.  This involves
      // a double loop to integrate the test functions (i) against
      // the trial functions (j). Note that the Jacobian depends
      // on the current solution x, which we access using the soln
      // vector.
      //
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          Gradient grad_u;
 
          for (unsigned int i=0; i<n_dofs; i++)
            grad_u += dphi[i][qp]*soln(dof_indices[i]);
 
          const Number
            sa = 1. + grad_u*grad_u,
            K  = 1. / std::sqrt(sa),
            dK = -K / sa;
 
          for (unsigned int i=0; i<n_dofs; i++)
            for (unsigned int j=0; j<n_dofs; j++)
              Ke(i,j) += JxW[qp]*(
                                  K * (dphi[i][qp]*dphi[j][qp]) +
                                  dK * (grad_u*dphi[j][qp]) * (grad_u*dphi[i][qp]) +
                                  _kappa * phi[i][qp] * phi[j][qp]
                                  );
        }
 
      dof_map.constrain_element_matrix (Ke, dof_indices);
 
      // Add the element matrix to the system Jacobian.
      jacobian.add_matrix (Ke, dof_indices);
    }
 
  // That's it.
}

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::FEGenericBase< OutputType >::build(), dim, libMesh::FIFTH, libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::DenseMatrix< T >::resize(), and std::sqrt().

jacobian() [2/2]

virtual void libMesh::NonlinearImplicitSystem::ComputeJacobian::jacobian ( const NumericVector< Number > &  X, SparseMatrix< Number > &  J, sys_type &  S  )

pure virtualinherited

Jacobian function.

This function will be called to compute the jacobian and must be implemented by the user in a derived class.

Referenced by libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), and libMesh::libmesh_petsc_snes_jacobian().

postcheck() [1/2]

void LaplaceYoung::postcheck ( const NumericVector< Number > &  old_soln, NumericVector< Number > &  search_direction, NumericVector< Number > &  new_soln, bool &  changed_search_direction, bool &  changed_new_soln, NonlinearImplicitSystem &  S  )

virtual

Function which performs a postcheck on the solution.

In this example, we take a "damped" Newton step defined by:

u_new = u_old + gamma*delta_u

where delta_u is the search direction, and 0 < gamma <= 1 is a damping parameter. gamma=1 corresponds to a full Newton step.

This is really for demonstration purposes only, as it just degrades the rate of nonlinear convergence in this particular example.

Definition at line 592 of file miscellaneous_ex3.C.

{
  // Back up along the search direction by some amount.  Since Newton
  // already works well for this problem, the only affect of this is
  // to degrade the rate of convergence.
  //
  // The minus sign is due to the sign of the "search_direction"
  // vector which comes in from the Newton solve. The RHS of the
  // nonlinear system, i.e. the residual, is generally not multiplied
  // by -1, so the solution vector, i.e. the search_direction, has a
  // factor -1.
  if (_gamma != 1.0)
    {
      new_soln = old_soln;
      new_soln.add(-_gamma, search_direction);
      changed_new_soln = true;
    }
  else
    changed_new_soln = false;
}

References libMesh::NumericVector< T >::add().

postcheck() [2/2]

virtual void libMesh::NonlinearImplicitSystem::ComputePostCheck::postcheck ( const NumericVector< Number > &  old_soln, NumericVector< Number > &  search_direction, NumericVector< Number > &  new_soln, bool &  changed_search_direction, bool &  changed_new_soln, sys_type &  S  )

pure virtualinherited

This interface, which is inspired by PETSc's, passes the user: .) A constant reference to the "old" solution vector (previous Newton iterate), .) A pointer to the search direction (update) vector, and .) The "proposed" new solution vector.

The user can then modify either the search direction or the new solution vector according to their own application-specific criteria. If they do, they must set the associated boolean references appropriately.

Referenced by libMesh::libmesh_petsc_snes_postcheck().

residual() [1/2]

void LaplaceYoung::residual ( const NumericVector< Number > &  X, NumericVector< Number > &  R, NonlinearImplicitSystem &  S  )

virtual

Function which computes the residual.

Definition at line 294 of file miscellaneous_ex3.C.

{
  EquationSystems & es = sys.get_equation_systems();
 
  // 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();
  libmesh_assert_equal_to (dim, 2);
 
  // Get a reference to the NonlinearImplicitSystem we are solving
  NonlinearImplicitSystem & system =
    es.get_system<NonlinearImplicitSystem>("Laplace-Young");
 
  // 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.  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.
  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
 
  // A 5th order Gauss quadrature rule for numerical integration.
  QGauss qrule (dim, FIFTH);
 
  // Tell the finite element object to use our quadrature rule.
  fe->attach_quadrature_rule (&qrule);
 
  // Declare a special finite element object for
  // boundary integration.
  std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
 
  // Boundary integration requires one quadrature rule,
  // with dimensionality one less than the dimensionality
  // of the element.
  QGauss qface(dim-1, FIFTH);
 
  // Tell the finite element object to use our
  // quadrature rule.
  fe_face->attach_quadrature_rule (&qface);
 
  // Here we define some references to cell-specific data that
  // will be used to assemble the linear system.
  // We begin with the element Jacobian * quadrature weight at each
  // integration point.
  const std::vector<Real> & JxW = fe->get_JxW();
 
  // The element shape functions evaluated at the quadrature points.
  const std::vector<std::vector<Real>> & phi = fe->get_phi();
 
  // The element shape function gradients evaluated at the quadrature
  // points.
  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
 
  // Define data structures to contain the residual contributions
  DenseVector<Number> Re;
 
  // 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;
 
  // Now we will loop over all the active elements in the mesh which
  // are local to this processor.
  // We will compute the element residual.
  residual.zero();
 
  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);
 
      // 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).
      const unsigned int n_dofs =
        cast_int<unsigned int>(dof_indices.size());
      Re.resize (n_dofs);
 
      // Now we will build the residual. This involves
      // the construction of the matrix K and multiplication of it
      // with the current solution x. We rearrange this into two loops:
      // In the first, we calculate only the contribution of
      // K_ij*x_j which is independent of the row i. In the second loops,
      // we multiply with the row-dependent part and add it to the element
      // residual.
 
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        {
          Number u = 0;
          Gradient grad_u;
 
          for (unsigned int j=0; j<n_dofs; j++)
            {
              u      += phi[j][qp]*soln(dof_indices[j]);
              grad_u += dphi[j][qp]*soln(dof_indices[j]);
            }
 
          const Number K = 1./std::sqrt(1. + grad_u*grad_u);
 
          for (unsigned int i=0; i<n_dofs; i++)
            Re(i) += JxW[qp]*(
                              K*(dphi[i][qp]*grad_u) +
                              _kappa*phi[i][qp]*u
                              );
        }
 
      // At this point the interior element integration has
      // been completed.  However, we have not yet addressed
      // boundary conditions.
 
      // The following loops 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<Real>> & phi_face = fe_face->get_phi();
 
            // The Jacobian * Quadrature Weight at the quadrature
            // points on the face.
            const std::vector<Real> & JxW_face = fe_face->get_JxW();
 
            // Compute the shape function values on the element face.
            fe_face->reinit(elem, side);
 
            // Loop over the face quadrature points for integration.
            for (unsigned int qp=0; qp<qface.n_points(); qp++)
              {
                // This is the right-hand-side contribution (f),
                // which has to be subtracted from the current residual
                for (unsigned int i=0; i<n_dofs; i++)
                  Re(i) -= JxW_face[qp]*_sigma*phi_face[i][qp];
              }
          }
 
      dof_map.constrain_element_vector (Re, dof_indices);
      residual.add_vector (Re, dof_indices);
    }
 
  // That's it.
}

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), dim, libMesh::FIFTH, libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::DenseVector< T >::resize(), std::sqrt(), and libMesh::NumericVector< T >::zero().

residual() [2/2]

virtual void libMesh::NonlinearImplicitSystem::ComputeResidual::residual ( const NumericVector< Number > &  X, NumericVector< Number > &  R, sys_type &  S  )

pure virtualinherited

Residual function.

This function will be called to compute the residual and must be implemented by the user in a derived class.

Referenced by libMesh::Problem_Interface::computeF(), libMesh::libmesh_petsc_snes_fd_residual(), libMesh::libmesh_petsc_snes_mffd_residual(), and libMesh::libmesh_petsc_snes_residual().

Member Data Documentation

_gamma

Real LaplaceYoung::_gamma

private

Definition at line 128 of file miscellaneous_ex3.C.

_kappa

Real LaplaceYoung::_kappa

private

Definition at line 124 of file miscellaneous_ex3.C.

_sigma

Real LaplaceYoung::_sigma

private

Definition at line 125 of file miscellaneous_ex3.C.

---

The documentation for this class was generated from the following file:

• examples/miscellaneous/miscellaneous_ex3/miscellaneous_ex3.C