libMesh
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...

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.

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.
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
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  {
559
560  for (std::size_t i=0; i<phi.size(); i++)
562
563  const Number
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]) +
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.
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.

 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.

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

 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.

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.
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
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
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;
405
406  for (std::size_t j=0; j<phi.size(); j++)
407  {
408  u += phi[j][qp]*soln(dof_indices[j]);
410  }
411
413
414  for (std::size_t i=0; i<phi.size(); i++)
415  Re(i) += JxW[qp]*(
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
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);
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: