libMesh
Public Member Functions | Private Attributes | List of all members
LaplaceYoung Class Referenceabstract

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

Inheritance diagram for LaplaceYoung:
[legend]

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

Definition at line 83 of file miscellaneous_ex3.C.

83  :
84  _kappa(1.),
85  _sigma(0.2),
86  _gamma(1.0)
87  {}

Member Function Documentation

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

Function which computes the jacobian.

Definition at line 462 of file miscellaneous_ex3.C.

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(), and libMesh::DenseMatrix< T >::resize().

465 {
466  // Get a reference to the equation system.
468 
469  // Get a constant reference to the mesh object.
470  const MeshBase & mesh = es.get_mesh();
471 
472  // The dimension that we are running
473  const unsigned int dim = mesh.mesh_dimension();
474 
475  // Get a reference to the NonlinearImplicitSystem we are solving
476  NonlinearImplicitSystem & system =
477  es.get_system<NonlinearImplicitSystem>("Laplace-Young");
478 
479  // A reference to the DofMap object for this system. The DofMap
480  // object handles the index translation from node and element numbers
481  // to degree of freedom numbers. We will talk more about the DofMap
482  // in future examples.
483  const DofMap & dof_map = system.get_dof_map();
484 
485  // Get a constant reference to the Finite Element type
486  // for the first (and only) variable in the system.
487  FEType fe_type = dof_map.variable_type(0);
488 
489  // Build a Finite Element object of the specified type. Since the
490  // FEBase::build() member dynamically creates memory we will
491  // store the object as a UniquePtr<FEBase>. This can be thought
492  // of as a pointer that will clean up after itself.
493  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
494 
495  // A 5th order Gauss quadrature rule for numerical integration.
496  QGauss qrule (dim, FIFTH);
497 
498  // Tell the finite element object to use our quadrature rule.
499  fe->attach_quadrature_rule (&qrule);
500 
501  // Here we define some references to cell-specific data that
502  // will be used to assemble the linear system.
503  // We begin with the element Jacobian * quadrature weight at each
504  // integration point.
505  const std::vector<Real> & JxW = fe->get_JxW();
506 
507  // The element shape functions evaluated at the quadrature points.
508  const std::vector<std::vector<Real>> & phi = fe->get_phi();
509 
510  // The element shape function gradients evaluated at the quadrature
511  // points.
512  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
513 
514  // Define data structures to contain the Jacobian element matrix.
515  // Following basic finite element terminology we will denote these
516  // "Ke". More detail is in example 3.
518 
519  // This vector will hold the degree of freedom indices for
520  // the element. These define where in the global system
521  // the element degrees of freedom get mapped.
522  std::vector<dof_id_type> dof_indices;
523 
524  // Now we will loop over all the active elements in the mesh which
525  // are local to this processor.
526  // We will compute the element Jacobian contribution.
527  for (const auto & elem : mesh.active_local_element_ptr_range())
528  {
529  // Get the degree of freedom indices for the
530  // current element. These define where in the global
531  // matrix and right-hand-side this element will
532  // contribute to.
533  dof_map.dof_indices (elem, dof_indices);
534 
535  // Compute the element-specific data for the current
536  // element. This involves computing the location of the
537  // quadrature points (q_point) and the shape functions
538  // (phi, dphi) for the current element.
539  fe->reinit (elem);
540 
541  // Zero the element Jacobian before
542  // summing them. We use the resize member here because
543  // the number of degrees of freedom might have changed from
544  // the last element. Note that this will be the case if the
545  // element type is different (i.e. the last element was a
546  // triangle, now we are on a quadrilateral).
547  Ke.resize (dof_indices.size(),
548  dof_indices.size());
549 
550  // Now we will build the element Jacobian. This involves
551  // a double loop to integrate the test functions (i) against
552  // the trial functions (j). Note that the Jacobian depends
553  // on the current solution x, which we access using the soln
554  // vector.
555  //
556  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
557  {
558  Gradient grad_u;
559 
560  for (std::size_t i=0; i<phi.size(); i++)
561  grad_u += dphi[i][qp]*soln(dof_indices[i]);
562 
563  const Number
564  sa = 1. + grad_u*grad_u,
565  K = 1. / std::sqrt(sa),
566  dK = -K / sa;
567 
568  for (std::size_t i=0; i<phi.size(); i++)
569  for (std::size_t j=0; j<phi.size(); j++)
570  Ke(i,j) += JxW[qp]*(
571  K * (dphi[i][qp]*dphi[j][qp]) +
572  dK * (grad_u*dphi[j][qp]) * (grad_u*dphi[i][qp]) +
573  _kappa * phi[i][qp] * phi[j][qp]
574  );
575  }
576 
577  dof_map.constrain_element_matrix (Ke, dof_indices);
578 
579  // Add the element matrix to the system Jacobian.
580  jacobian.add_matrix (Ke, dof_indices);
581  }
582 
583  // That's it.
584 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
unsigned int dim
ImplicitSystem & sys
MeshBase & mesh
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
const DofMap & get_dof_map() const
Definition: system.h:2030
This class provides a specific system class.
virtual void jacobian(const NumericVector< Number > &X, SparseMatrix< Number > &J, NonlinearImplicitSystem &S)
Function which computes the jacobian.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
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::__libmesh_petsc_snes_jacobian(), libMesh::Problem_Interface::computeJacobian(), and libMesh::Problem_Interface::computePreconditioner().

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 589 of file miscellaneous_ex3.C.

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

595 {
596  // Back up along the search direction by some amount. Since Newton
597  // already works well for this problem, the only affect of this is
598  // to degrade the rate of convergence.
599  //
600  // The minus sign is due to the sign of the "search_direction"
601  // vector which comes in from the Newton solve. The RHS of the
602  // nonlinear system, i.e. the residual, is generally not multiplied
603  // by -1, so the solution vector, i.e. the search_direction, has a
604  // factor -1.
605  if (_gamma != 1.0)
606  {
607  new_soln = old_soln;
608  new_soln.add(-_gamma, search_direction);
609  changed_new_soln = true;
610  }
611  else
612  changed_new_soln = false;
613 }
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
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().

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::__libmesh_petsc_snes_residual(), libMesh::Problem_Interface::computeF(), and libMesh::NonlinearImplicitSystem::ComputeResidual::~ComputeResidual().

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.

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(), libmesh_nullptr, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::QBase::n_points(), libMesh::DenseVector< T >::resize(), side, and libMesh::NumericVector< T >::zero().

297 {
299 
300  // Get a constant reference to the mesh object.
301  const MeshBase & mesh = es.get_mesh();
302 
303  // The dimension that we are running
304  const unsigned int dim = mesh.mesh_dimension();
305  libmesh_assert_equal_to (dim, 2);
306 
307  // Get a reference to the NonlinearImplicitSystem we are solving
308  NonlinearImplicitSystem & system =
309  es.get_system<NonlinearImplicitSystem>("Laplace-Young");
310 
311  // A reference to the DofMap object for this system. The DofMap
312  // object handles the index translation from node and element numbers
313  // to degree of freedom numbers. We will talk more about the DofMap
314  // in future examples.
315  const DofMap & dof_map = system.get_dof_map();
316 
317  // Get a constant reference to the Finite Element type
318  // for the first (and only) variable in the system.
319  FEType fe_type = dof_map.variable_type(0);
320 
321  // Build a Finite Element object of the specified type. Since the
322  // FEBase::build() member dynamically creates memory we will
323  // store the object as a UniquePtr<FEBase>. This can be thought
324  // of as a pointer that will clean up after itself.
325  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
326 
327  // A 5th order Gauss quadrature rule for numerical integration.
328  QGauss qrule (dim, FIFTH);
329 
330  // Tell the finite element object to use our quadrature rule.
331  fe->attach_quadrature_rule (&qrule);
332 
333  // Declare a special finite element object for
334  // boundary integration.
335  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
336 
337  // Boundary integration requires one quadrature rule,
338  // with dimensionality one less than the dimensionality
339  // of the element.
340  QGauss qface(dim-1, FIFTH);
341 
342  // Tell the finite element object to use our
343  // quadrature rule.
344  fe_face->attach_quadrature_rule (&qface);
345 
346  // Here we define some references to cell-specific data that
347  // will be used to assemble the linear system.
348  // We begin with the element Jacobian * quadrature weight at each
349  // integration point.
350  const std::vector<Real> & JxW = fe->get_JxW();
351 
352  // The element shape functions evaluated at the quadrature points.
353  const std::vector<std::vector<Real>> & phi = fe->get_phi();
354 
355  // The element shape function gradients evaluated at the quadrature
356  // points.
357  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
358 
359  // Define data structures to contain the residual contributions
361 
362  // This vector will hold the degree of freedom indices for
363  // the element. These define where in the global system
364  // the element degrees of freedom get mapped.
365  std::vector<dof_id_type> dof_indices;
366 
367  // Now we will loop over all the active elements in the mesh which
368  // are local to this processor.
369  // We will compute the element residual.
370  residual.zero();
371 
372  for (const auto & elem : mesh.active_local_element_ptr_range())
373  {
374  // Get the degree of freedom indices for the
375  // current element. These define where in the global
376  // matrix and right-hand-side this element will
377  // contribute to.
378  dof_map.dof_indices (elem, dof_indices);
379 
380  // Compute the element-specific data for the current
381  // element. This involves computing the location of the
382  // quadrature points (q_point) and the shape functions
383  // (phi, dphi) for the current element.
384  fe->reinit (elem);
385 
386  // We use the resize member here because
387  // the number of degrees of freedom might have changed from
388  // the last element. Note that this will be the case if the
389  // element type is different (i.e. the last element was a
390  // triangle, now we are on a quadrilateral).
391  Re.resize (dof_indices.size());
392 
393  // Now we will build the residual. This involves
394  // the construction of the matrix K and multiplication of it
395  // with the current solution x. We rearrange this into two loops:
396  // In the first, we calculate only the contribution of
397  // K_ij*x_j which is independent of the row i. In the second loops,
398  // we multiply with the row-dependent part and add it to the element
399  // residual.
400 
401  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
402  {
403  Number u = 0;
404  Gradient grad_u;
405 
406  for (std::size_t j=0; j<phi.size(); j++)
407  {
408  u += phi[j][qp]*soln(dof_indices[j]);
409  grad_u += dphi[j][qp]*soln(dof_indices[j]);
410  }
411 
412  const Number K = 1./std::sqrt(1. + grad_u*grad_u);
413 
414  for (std::size_t i=0; i<phi.size(); i++)
415  Re(i) += JxW[qp]*(
416  K*(dphi[i][qp]*grad_u) +
417  _kappa*phi[i][qp]*u
418  );
419  }
420 
421  // At this point the interior element integration has
422  // been completed. However, we have not yet addressed
423  // boundary conditions.
424 
425  // The following loops over the sides of the element.
426  // If the element has no neighbor on a side then that
427  // side MUST live on a boundary of the domain.
428  for (auto side : elem->side_index_range())
429  if (elem->neighbor_ptr(side) == libmesh_nullptr)
430  {
431  // The value of the shape functions at the quadrature
432  // points.
433  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
434 
435  // The Jacobian * Quadrature Weight at the quadrature
436  // points on the face.
437  const std::vector<Real> & JxW_face = fe_face->get_JxW();
438 
439  // Compute the shape function values on the element face.
440  fe_face->reinit(elem, side);
441 
442  // Loop over the face quadrature points for integration.
443  for (unsigned int qp=0; qp<qface.n_points(); qp++)
444  {
445  // This is the right-hand-side contribution (f),
446  // which has to be subtracted from the current residual
447  for (std::size_t i=0; i<phi_face.size(); i++)
448  Re(i) -= JxW_face[qp]*_sigma*phi_face[i][qp];
449  }
450  }
451 
452  dof_map.constrain_element_vector (Re, dof_indices);
453  residual.add_vector (Re, dof_indices);
454  }
455 
456  // That's it.
457 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
unsigned int dim
ImplicitSystem & sys
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
const DofMap & get_dof_map() const
Definition: system.h:2030
This class provides a specific system class.
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
This class implements specific orders of Gauss quadrature.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
virtual void residual(const NumericVector< Number > &X, NumericVector< Number > &R, NonlinearImplicitSystem &S)
Function which computes the residual.

Member Data Documentation

Real LaplaceYoung::_gamma
private

Definition at line 128 of file miscellaneous_ex3.C.

Real LaplaceYoung::_kappa
private

Definition at line 124 of file miscellaneous_ex3.C.

Real LaplaceYoung::_sigma
private

Definition at line 125 of file miscellaneous_ex3.C.


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