libMesh
Public Member Functions | Private Attributes | List of all members
LinearElasticityWithContact Class Reference

This class encapsulate all functionality required for assembling and solving a linear elastic model with contact. More...

#include <linear_elasticity_with_contact.h>

Inheritance diagram for LinearElasticityWithContact:
[legend]

Public Member Functions

 LinearElasticityWithContact (NonlinearImplicitSystem &sys_in, Real contact_penalty_in)
 Constructor. More...
 
void set_contact_penalty (Real contact_penalty_in)
 Update the penalty parameter. More...
 
Real get_contact_penalty () const
 Get the penalty parameter. More...
 
Real kronecker_delta (unsigned int i, unsigned int j)
 Kronecker delta function. More...
 
Real elasticity_tensor (Real young_modulus, Real poisson_ratio, unsigned int i, unsigned int j, unsigned int k, unsigned int l)
 Evaluate the fourth order tensor (C_ijkl) that relates stress to strain. More...
 
void move_mesh (MeshBase &input_mesh, const NumericVector< Number > &input_solution)
 Move the mesh nodes of input_mesh based on the displacement field in input_solution. More...
 
void initialize_contact_load_paths ()
 Set up the load paths on the contact surfaces. More...
 
void add_contact_edge_elements ()
 Add edge elements into the mesh based on the contact load paths. More...
 
virtual void residual_and_jacobian (const NumericVector< Number > &soln, NumericVector< Number > *residual, SparseMatrix< Number > *jacobian, NonlinearImplicitSystem &)
 Evaluate the Jacobian of the nonlinear system. More...
 
void compute_stresses ()
 Compute the Cauchy stress for the current solution. More...
 
std::pair< Real, Real > update_lambdas ()
 Update the lambda parameters in the augmented Lagrangian method. More...
 
std::pair< Real, Real > get_least_and_max_gap_function ()
 

Private Attributes

NonlinearImplicitSystem_sys
 Keep a reference to the NonlinearImplicitSystem. More...
 
Real _contact_penalty
 Penalize overlapping elements. More...
 
std::map< dof_id_type, Real > _lambda_plus_penalty_values
 Store the intermediate values of lambda plus penalty. More...
 
std::map< dof_id_type, Real > _lambdas
 Augmented Lagrangian values at each contact node. More...
 
std::map< dof_id_type, dof_id_type > _contact_node_map
 This provides a map between contact nodes. More...
 

Detailed Description

This class encapsulate all functionality required for assembling and solving a linear elastic model with contact.

Definition at line 31 of file linear_elasticity_with_contact.h.

Constructor & Destructor Documentation

LinearElasticityWithContact::LinearElasticityWithContact ( NonlinearImplicitSystem sys_in,
Real  contact_penalty_in 
)

Constructor.

Definition at line 33 of file linear_elasticity_with_contact.C.

34  :
35  _sys(sys_in),
36  _contact_penalty(contact_penalty_in)
37 {
38 }
Real _contact_penalty
Penalize overlapping elements.
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.

Member Function Documentation

void LinearElasticityWithContact::add_contact_edge_elements ( )

Add edge elements into the mesh based on the contact load paths.

This ensure proper parallel communication of data, e.g. if a node on one side of the contact surface has a constraint on it, then adding contact elements into the mesh ensures that the constraint will be enforced properly.

Definition at line 207 of file linear_elasticity_with_contact.C.

References _contact_node_map, _sys, libMesh::MeshBase::add_elem(), libMesh::System::get_mesh(), mesh, libMesh::MeshBase::node_ref(), libMesh::MeshBase::prepare_for_use(), libMesh::Elem::set_node(), and libMesh::Elem::subdomain_id().

208 {
209  MeshBase & mesh = _sys.get_mesh();
210 
211  std::map<dof_id_type, dof_id_type>::iterator it = _contact_node_map.begin();
212  std::map<dof_id_type, dof_id_type>::iterator it_end = _contact_node_map.end();
213  for( ; it != it_end ; ++it )
214  {
215  dof_id_type master_node_id = it->first;
216  dof_id_type slave_node_id = it->second;
217 
218  Node & master_node = mesh.node_ref(master_node_id);
219  Node & slave_node = mesh.node_ref(slave_node_id);
220 
221  Elem * connector_elem = mesh.add_elem (new Edge2);
222  connector_elem->set_node(0) = &master_node;
223  connector_elem->set_node(1) = &slave_node;
224 
225  connector_elem->subdomain_id() = 10;
226  }
227 
228  mesh.prepare_for_use();
229 }
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1941
A Node is like a Point, but with more information.
Definition: node.h:52
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:420
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
This is the MeshBase class.
Definition: mesh_base.h:68
const MeshBase & get_mesh() const
Definition: system.h:2014
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:174
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
The Edge2 is an element in 1D composed of 2 nodes.
Definition: edge_edge2.h:43
uint8_t dof_id_type
Definition: id_types.h:64
void LinearElasticityWithContact::compute_stresses ( )

Compute the Cauchy stress for the current solution.

Definition at line 481 of file linear_elasticity_with_contact.C.

References _sys, libMesh::DenseMatrix< T >::add(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::EDGE2, elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::System::n_vars(), libMesh::EquationSystems::parameters, std::pow(), libMesh::Real, libMesh::DenseMatrix< T >::resize(), libMesh::DenseMatrix< T >::scale(), libMesh::System::solution, libMesh::System::update(), libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

482 {
484  const Real young_modulus = es.parameters.get<Real>("young_modulus");
485  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
486 
487  const MeshBase & mesh = _sys.get_mesh();
488  const unsigned int dim = mesh.mesh_dimension();
489 
490  unsigned int displacement_vars[3];
491  displacement_vars[0] = _sys.variable_number ("u");
492  displacement_vars[1] = _sys.variable_number ("v");
493  displacement_vars[2] = _sys.variable_number ("w");
494  const unsigned int u_var = _sys.variable_number ("u");
495 
496  const DofMap & dof_map = _sys.get_dof_map();
497  FEType fe_type = dof_map.variable_type(u_var);
498  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
499  QGauss qrule (dim, fe_type.default_quadrature_order());
500  fe->attach_quadrature_rule (&qrule);
501 
502  const std::vector<Real> & JxW = fe->get_JxW();
503  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
504 
505  // Also, get a reference to the ExplicitSystem
506  ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
507  const DofMap & stress_dof_map = stress_system.get_dof_map();
508  unsigned int sigma_vars[6];
509  sigma_vars[0] = stress_system.variable_number ("sigma_00");
510  sigma_vars[1] = stress_system.variable_number ("sigma_01");
511  sigma_vars[2] = stress_system.variable_number ("sigma_02");
512  sigma_vars[3] = stress_system.variable_number ("sigma_11");
513  sigma_vars[4] = stress_system.variable_number ("sigma_12");
514  sigma_vars[5] = stress_system.variable_number ("sigma_22");
515  unsigned int vonMises_var = stress_system.variable_number ("vonMises");
516 
517  // Storage for the stress dof indices on each element
518  std::vector<std::vector<dof_id_type>> dof_indices_var(_sys.n_vars());
519  std::vector<dof_id_type> stress_dof_indices_var;
520 
521  // To store the stress tensor on each element
522  DenseMatrix<Number> elem_avg_stress_tensor(3, 3);
523 
524  for (const auto & elem : mesh.active_local_element_ptr_range())
525  {
526  if( elem->type() == EDGE2 )
527  {
528  // We do not compute stress on the contact connector elements.
529  continue;
530  }
531 
532  for (unsigned int var=0; var<3; var++)
533  dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
534 
535  const unsigned int n_var_dofs = dof_indices_var[0].size();
536 
537  fe->reinit (elem);
538 
539  // clear the stress tensor
540  elem_avg_stress_tensor.resize(3, 3);
541 
542  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
543  {
544  // Row is variable u1, u2, or u3, column is x, y, or z
545  DenseMatrix<Number> grad_u(3, 3);
546  for (unsigned int var_i=0; var_i<3; var_i++)
547  for (unsigned int var_j=0; var_j<3; var_j++)
548  for (unsigned int j=0; j<n_var_dofs; j++)
549  grad_u(var_i, var_j) +=
550  dphi[j][qp](var_j) * _sys.current_solution(dof_indices_var[var_i][j]);
551 
552  DenseMatrix<Number> stress_tensor(3, 3);
553  for (unsigned int i=0; i<3; i++)
554  for (unsigned int j=0; j<3; j++)
555  for (unsigned int k=0; k<3; k++)
556  for (unsigned int l=0; l<3; l++)
557  stress_tensor(i,j) +=
558  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * grad_u(k,l);
559 
560  // We want to plot the average stress on each element, hence
561  // we integrate stress_tensor
562  elem_avg_stress_tensor.add(JxW[qp], stress_tensor);
563  }
564 
565  // Get the average stress per element by dividing by volume
566  elem_avg_stress_tensor.scale(1./elem->volume());
567 
568  // load elem_sigma data into stress_system
569  unsigned int stress_var_index = 0;
570  for (unsigned int i=0; i<3; i++)
571  for (unsigned int j=i; j<3; j++)
572  {
573  stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[stress_var_index]);
574 
575  // We are using CONSTANT MONOMIAL basis functions, hence we only need to get
576  // one dof index per variable
577  dof_id_type dof_index = stress_dof_indices_var[0];
578 
579  if ((stress_system.solution->first_local_index() <= dof_index) &&
580  (dof_index < stress_system.solution->last_local_index()))
581  stress_system.solution->set(dof_index, elem_avg_stress_tensor(i,j));
582 
583  stress_var_index++;
584  }
585 
586  // Also, the von Mises stress
587  Number vonMises_value = std::sqrt(0.5*(pow(elem_avg_stress_tensor(0,0) - elem_avg_stress_tensor(1,1), 2.) +
588  pow(elem_avg_stress_tensor(1,1) - elem_avg_stress_tensor(2,2), 2.) +
589  pow(elem_avg_stress_tensor(2,2) - elem_avg_stress_tensor(0,0), 2.) +
590  6.*(pow(elem_avg_stress_tensor(0,1), 2.) +
591  pow(elem_avg_stress_tensor(1,2), 2.) +
592  pow(elem_avg_stress_tensor(2,0), 2.))));
593 
594  stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
595  dof_id_type dof_index = stress_dof_indices_var[0];
596 
597  if ((stress_system.solution->first_local_index() <= dof_index) &&
598  (dof_index < stress_system.solution->last_local_index()))
599  stress_system.solution->set(dof_index, vonMises_value);
600  }
601 
602  // Should call close and update when we set vector entries directly
603  stress_system.solution->close();
604  stress_system.update();
605 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
MeshBase & mesh
const T & get(const std::string &) const
Definition: parameters.h:430
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
Order default_quadrature_order() const
Definition: fe_type.h:332
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
double pow(double a, int b)
Real elasticity_tensor(Real young_modulus, Real poisson_ratio, unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
This class implements specific orders of Gauss quadrature.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
Real LinearElasticityWithContact::elasticity_tensor ( Real  young_modulus,
Real  poisson_ratio,
unsigned int  i,
unsigned int  j,
unsigned int  k,
unsigned int  l 
)

Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.

Definition at line 56 of file linear_elasticity_with_contact.C.

References kronecker_delta(), and libMesh::Real.

Referenced by compute_stresses(), and residual_and_jacobian().

62 {
63  // Define the Lame constants
64  const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
65  const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
66 
67  return lambda_1 * kronecker_delta(i, j) * kronecker_delta(k, l) +
68  lambda_2 * (kronecker_delta(i, k) * kronecker_delta(j, l) + kronecker_delta(i, l) * kronecker_delta(j, k));
69 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real kronecker_delta(unsigned int i, unsigned int j)
Kronecker delta function.
Real LinearElasticityWithContact::get_contact_penalty ( ) const

Get the penalty parameter.

Definition at line 45 of file linear_elasticity_with_contact.C.

References _contact_penalty.

46 {
47  return _contact_penalty;
48 }
Real _contact_penalty
Penalize overlapping elements.
std::pair< Real, Real > LinearElasticityWithContact::get_least_and_max_gap_function ( )
Returns
the least and max gap function values for the current solution.

Definition at line 638 of file linear_elasticity_with_contact.C.

References _contact_node_map, _sys, libMesh::MeshBase::clone(), libMesh::System::get_mesh(), std::max(), std::min(), move_mesh(), libMesh::Real, and libMesh::System::solution.

639 {
640  UniquePtr<MeshBase> mesh_clone = _sys.get_mesh().clone();
641  move_mesh(*mesh_clone, *_sys.solution);
642 
643  Real least_value = std::numeric_limits<Real>::max();
644  Real max_value = std::numeric_limits<Real>::min();
645 
646  std::map<dof_id_type, dof_id_type>::iterator it = _contact_node_map.begin();
647  std::map<dof_id_type, dof_id_type>::iterator it_end = _contact_node_map.end();
648  for ( ; it != it_end; ++it)
649  {
650  dof_id_type lower_point_id = it->first;
651  dof_id_type upper_point_id = it->second;
652 
653  Point upper_to_lower;
654  {
655  Point lower_node_moved = mesh_clone->point(lower_point_id);
656  Point upper_node_moved = mesh_clone->point(upper_point_id);
657 
658  upper_to_lower = (lower_node_moved - upper_node_moved);
659  }
660 
661  // Set the contact force direction. Usually this would be calculated
662  // separately on each master node, but here we just set it to (0, 0, 1)
663  // everywhere.
664  Point contact_force_direction(0., 0., 1.);
665 
666  // gap_function > 0. means that contact has been detected
667  // gap_function < 0. means that we have a gap
668  // This sign convention matches Simo & Laursen (1992).
669  Real gap_function = upper_to_lower * contact_force_direction;
670 
671  if (gap_function < least_value)
672  least_value = gap_function;
673 
674  if (gap_function > max_value)
675  max_value = gap_function;
676  }
677 
678  return std::make_pair(least_value, max_value);
679 }
long double max(long double a, double b)
void move_mesh(MeshBase &input_mesh, const NumericVector< Number > &input_solution)
Move the mesh nodes of input_mesh based on the displacement field in input_solution.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
const MeshBase & get_mesh() const
Definition: system.h:2014
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
virtual UniquePtr< MeshBase > clone() const =0
Virtual "copy constructor".
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
uint8_t dof_id_type
Definition: id_types.h:64
void LinearElasticityWithContact::initialize_contact_load_paths ( )

Set up the load paths on the contact surfaces.

Definition at line 142 of file linear_elasticity_with_contact.C.

References _contact_node_map, _lambdas, _sys, libMesh::MeshBase::active_element_ptr_range(), distance(), libMesh::MeshBase::get_boundary_info(), libMesh::System::get_mesh(), libMesh::BoundaryInfo::has_boundary_id(), libMesh::libmesh_assert(), libmesh_nullptr, std::max(), mesh, libMesh::MeshBase::point(), libMesh::Real, and side.

143 {
144  const MeshBase & mesh = _sys.get_mesh();
145 
146  std::vector<dof_id_type> nodes_on_lower_surface;
147  std::vector<dof_id_type> nodes_on_upper_surface;
148 
149  _lambdas.clear();
150  for (const auto & elem : mesh.active_element_ptr_range())
151  for (auto side : elem->side_index_range())
152  if (elem->neighbor_ptr(side) == libmesh_nullptr)
153  {
154  bool on_lower_contact_surface =
155  mesh.get_boundary_info().has_boundary_id (elem, side, CONTACT_BOUNDARY_LOWER);
156 
157  bool on_upper_contact_surface =
158  mesh.get_boundary_info().has_boundary_id (elem, side, CONTACT_BOUNDARY_UPPER);
159 
160  if (on_lower_contact_surface && on_upper_contact_surface)
161  libmesh_error_msg("Should not be on both surfaces at the same time");
162 
163  if (on_lower_contact_surface || on_upper_contact_surface)
164  {
165  for (auto node_index : elem->node_index_range())
166  if (elem->is_node_on_side(node_index, side))
167  {
168  if (on_lower_contact_surface)
169  nodes_on_lower_surface.push_back(elem->node_id(node_index));
170  else
171  {
172  _lambdas[elem->node_id(node_index)] = 0.;
173  nodes_on_upper_surface.push_back(elem->node_id(node_index));
174  }
175  }
176  }
177  } // end if neighbor(side_) != libmesh_nullptr
178 
179  // In this example, we expect the number of upper and lower nodes to match
180  libmesh_assert(nodes_on_lower_surface.size() == nodes_on_upper_surface.size());
181 
182  // Do an N^2 search to match the contact nodes
183  _contact_node_map.clear();
184  for (std::size_t i=0; i<nodes_on_lower_surface.size(); i++)
185  {
186  dof_id_type lower_node_id = nodes_on_lower_surface[i];
187  Point p_lower = mesh.point(lower_node_id);
188 
189  Real min_distance = std::numeric_limits<Real>::max();
190 
191  for (std::size_t j=0; j<nodes_on_upper_surface.size(); j++)
192  {
193  dof_id_type upper_node_id = nodes_on_upper_surface[j];
194  Point p_upper = mesh.point(upper_node_id);
195 
196  Real distance = (p_upper - p_lower).norm();
197 
198  if (distance < min_distance)
199  {
200  _contact_node_map[lower_node_id] = upper_node_id;
201  min_distance = distance;
202  }
203  }
204  }
205 }
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
virtual const Point & point(const dof_id_type i) const =0
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
long double max(long double a, double b)
Real distance(const Point &p)
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
const MeshBase & get_mesh() const
Definition: system.h:2014
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
std::map< dof_id_type, Real > _lambdas
Augmented Lagrangian values at each contact node.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
uint8_t dof_id_type
Definition: id_types.h:64
Real LinearElasticityWithContact::kronecker_delta ( unsigned int  i,
unsigned int  j 
)

Kronecker delta function.

Definition at line 50 of file linear_elasticity_with_contact.C.

Referenced by elasticity_tensor().

52 {
53  return i == j ? 1. : 0.;
54 }
void LinearElasticityWithContact::move_mesh ( MeshBase input_mesh,
const NumericVector< Number > &  input_solution 
)

Move the mesh nodes of input_mesh based on the displacement field in input_solution.

Definition at line 71 of file linear_elasticity_with_contact.C.

References _sys, libMesh::MeshBase::active_element_ptr_range(), libMesh::ParallelObject::comm(), data, libMesh::DofMap::dof_indices(), libMesh::MeshBase::elem_ptr(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::NumericVector< T >::init(), libMesh::NumericVector< T >::localize(), libMesh::SERIAL, libMesh::NumericVector< T >::size(), value, libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

Referenced by get_least_and_max_gap_function(), and residual_and_jacobian().

73 {
74  // Maintain a set of node ids that we've encountered.
75  std::unordered_set<dof_id_type> encountered_node_ids;
76 
77  // Localize input_solution so that we have the data to move all
78  // elements (not just elements local to this processor).
79  UniquePtr<NumericVector<Number>> localized_input_solution =
80  NumericVector<Number>::build(input_solution.comm());
81 
82  localized_input_solution->init (input_solution.size(), false, SERIAL);
83  input_solution.localize(*localized_input_solution);
84 
85  for (const auto & elem : input_mesh.active_element_ptr_range())
86  {
87  Elem * orig_elem = _sys.get_mesh().elem_ptr(elem->id());
88 
89  for (auto node_id : elem->node_index_range())
90  {
91  Node & node = elem->node_ref(node_id);
92 
93  if (encountered_node_ids.find(node.id()) != encountered_node_ids.end())
94  continue;
95 
96  encountered_node_ids.insert(node.id());
97 
98  std::vector<std::string> uvw_names(3);
99  uvw_names[0] = "u";
100  uvw_names[1] = "v";
101  uvw_names[2] = "w";
102 
103  {
104  const Point master_point = elem->master_point(node_id);
105 
106  Point uvw;
107  for (std::size_t index=0; index<uvw_names.size(); index++)
108  {
109  const unsigned int var = _sys.variable_number(uvw_names[index]);
110  const FEType & fe_type = _sys.get_dof_map().variable_type(var);
111 
112  FEComputeData data (_sys.get_equation_systems(), master_point);
113 
114  FEInterface::compute_data(elem->dim(),
115  fe_type,
116  elem,
117  data);
118 
119  std::vector<dof_id_type> dof_indices_var;
120  _sys.get_dof_map().dof_indices (orig_elem, dof_indices_var, var);
121 
122  for (std::size_t i=0; i<dof_indices_var.size(); i++)
123  {
124  Number value = (*localized_input_solution)(dof_indices_var[i]) * data.shape[i];
125 
126 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
127  // We explicitly store the real part in uvw
128  uvw(index) += value.real();
129 #else
130  uvw(index) += value;
131 #endif
132  }
133  }
134 
135  // Update the node's location
136  node += uvw;
137  }
138  }
139  }
140 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
A Node is like a Point, but with more information.
Definition: node.h:52
virtual numeric_index_type size() const =0
virtual void init(const numeric_index_type, const numeric_index_type, const bool=false, const ParallelType=AUTOMATIC)=0
Change the dimension of the vector to N.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
class FEComputeData hides arbitrary data to be passed to and from children of FEBase through the FEIn...
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
Numeric vector.
Definition: dof_map.h:66
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
const Parallel::Communicator & comm() const
static const bool value
Definition: xdr_io.C:108
const EquationSystems & get_equation_systems() const
Definition: system.h:712
IterBase * data
Ideally this private member data should have protected access.
dof_id_type id() const
Definition: dof_object.h:632
virtual const Elem * elem_ptr(const dof_id_type i) const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
void LinearElasticityWithContact::residual_and_jacobian ( const NumericVector< Number > &  soln,
NumericVector< Number > *  residual,
SparseMatrix< Number > *  jacobian,
NonlinearImplicitSystem  
)
virtual

Evaluate the Jacobian of the nonlinear system.

Definition at line 231 of file linear_elasticity_with_contact.C.

References _contact_node_map, _contact_penalty, _lambda_plus_penalty_values, _lambdas, _sys, libMesh::SparseMatrix< T >::add(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::ParallelObject::comm(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), libMesh::DofObject::dof_number(), libMesh::EDGE2, elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::System::get_equation_systems(), libMesh::System::get_mesh(), mesh, move_mesh(), libMesh::System::number(), libMesh::EquationSystems::parameters, libMesh::Parallel::Communicator::rank(), libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::System::variable_number(), libMesh::DofMap::variable_type(), libMesh::SparseMatrix< T >::zero(), and libMesh::NumericVector< T >::zero().

235 {
237  const Real young_modulus = es.parameters.get<Real>("young_modulus");
238  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
239 
240  const MeshBase & mesh = _sys.get_mesh();
241  const unsigned int dim = mesh.mesh_dimension();
242 
243  const unsigned int u_var = _sys.variable_number ("u");
244 
245  DofMap & dof_map = _sys.get_dof_map();
246 
247  FEType fe_type = dof_map.variable_type(u_var);
248  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
249  QGauss qrule (dim, fe_type.default_quadrature_order());
250  fe->attach_quadrature_rule (&qrule);
251 
252  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
253  QGauss qface (dim-1, fe_type.default_quadrature_order());
254  fe_face->attach_quadrature_rule (&qface);
255 
256  UniquePtr<FEBase> fe_neighbor_face (FEBase::build(dim, fe_type));
257  fe_neighbor_face->attach_quadrature_rule (&qface);
258 
259  const std::vector<Real> & JxW = fe->get_JxW();
260  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
261 
262  if (jacobian)
263  jacobian->zero();
264 
265  if (residual)
266  residual->zero();
267 
268  // Do jacobian and residual assembly, including contact forces
270 
271  DenseSubVector<Number> Re_var[3] =
275 
277  DenseSubMatrix<Number> Ke_var[3][3] =
278  {
282  };
283 
284  std::vector<dof_id_type> dof_indices;
285  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
286 
287  for (const auto & elem : mesh.active_local_element_ptr_range())
288  {
289  if( elem->type() == EDGE2 )
290  {
291  // We do not do any assembly on the contact connector elements.
292  // The contact force assembly is handled in a separate loop.
293  continue;
294  }
295 
296  dof_map.dof_indices (elem, dof_indices);
297  for (unsigned int var=0; var<3; var++)
298  dof_map.dof_indices (elem, dof_indices_var[var], var);
299 
300  const unsigned int n_dofs = dof_indices.size();
301  const unsigned int n_var_dofs = dof_indices_var[0].size();
302 
303  fe->reinit (elem);
304 
305  Re.resize (n_dofs);
306  for (unsigned int var=0; var<3; var++)
307  Re_var[var].reposition (var*n_var_dofs, n_var_dofs);
308 
309  Ke.resize (n_dofs, n_dofs);
310  for (unsigned int var_i=0; var_i<3; var_i++)
311  for (unsigned int var_j=0; var_j<3; var_j++)
312  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
313 
314  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
315  {
316  // Row is variable u, v, or w column is x, y, or z
317  DenseMatrix<Number> grad_u(3, 3);
318  for (unsigned int var_i=0; var_i<3; var_i++)
319  for (unsigned int var_j=0; var_j<3; var_j++)
320  for (unsigned int j=0; j<n_var_dofs; j++)
321  grad_u(var_i, var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
322 
323  // - C_ijkl u_k,l v_i,j
324  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
325  for (unsigned int i=0; i<3; i++)
326  for (unsigned int j=0; j<3; j++)
327  for (unsigned int k=0; k<3; k++)
328  for (unsigned int l=0; l<3; l++)
329  Re_var[i](dof_i) -= JxW[qp] *
330  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * grad_u(k,l) * dphi[dof_i][qp](j);
331 
332  // assemble \int_Omega C_ijkl u_k,l v_i,j \dx
333  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
334  for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
335  for (unsigned int i=0; i<3; i++)
336  for (unsigned int j=0; j<3; j++)
337  for (unsigned int k=0; k<3; k++)
338  for (unsigned int l=0; l<3; l++)
339  Ke_var[i][k](dof_i, dof_j) -= JxW[qp] *
340  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * dphi[dof_j][qp](l) * dphi[dof_i][qp](j);
341  }
342 
343  dof_map.constrain_element_matrix_and_vector (Ke, Re, dof_indices);
344 
345  if (jacobian)
346  jacobian->add_matrix (Ke, dof_indices);
347 
348  if (residual)
349  residual->add_vector (Re, dof_indices);
350  }
351 
352  // Move a copy of the mesh based on the solution. This is used to get
353  // the contact forces. We could compute the contact forces based on the
354  // current displacement solution, which would not require a clone of the
355  // mesh. Avoiding the mesh clone would be important for production-scale
356  // contact solves, but for the sake of this example, using the clone is
357  // simple and fast enough.
358  UniquePtr<MeshBase> mesh_clone = mesh.clone();
359  move_mesh(*mesh_clone, soln);
360 
361  // Add contributions due to contact penalty forces. Only need to do this on
362  // one processor.
364 
365  std::map<dof_id_type, dof_id_type>::iterator it = _contact_node_map.begin();
366  std::map<dof_id_type, dof_id_type>::iterator it_end = _contact_node_map.end();
367 
368  for ( ; it != it_end; ++it)
369  {
370  dof_id_type lower_point_id = it->first;
371  dof_id_type upper_point_id = it->second;
372 
373  Point upper_to_lower;
374  {
375  Point lower_node_moved = mesh_clone->point(lower_point_id);
376  Point upper_node_moved = mesh_clone->point(upper_point_id);
377 
378  upper_to_lower = (lower_node_moved - upper_node_moved);
379  }
380 
381  // Set the contact force direction. Usually this would be calculated
382  // separately on each master node, but here we just set it to (0, 0, 1)
383  // everywhere.
384  Point contact_force_direction(0., 0., 1.);
385 
386  // gap_function > 0. means that contact has been detected
387  // gap_function < 0. means that we have a gap
388  // This sign convention matches Simo & Laursen (1992).
389  Real gap_function = upper_to_lower * contact_force_direction;
390 
391  // We use the sign of lambda_plus_penalty to determine whether or
392  // not we need to impose a contact force.
393  Real lambda_plus_penalty =
394  (_lambdas[upper_point_id] + gap_function * _contact_penalty);
395 
396  if (lambda_plus_penalty < 0.)
397  lambda_plus_penalty = 0.;
398 
399  // Store lambda_plus_penalty, we'll need to use it later to update _lambdas
400  _lambda_plus_penalty_values[upper_point_id] = lambda_plus_penalty;
401 
402  const Node & lower_node = mesh.node_ref(lower_point_id);
403  const Node & upper_node = mesh.node_ref(upper_point_id);
404 
405  std::vector<dof_id_type> dof_indices_on_lower_node(3);
406  std::vector<dof_id_type> dof_indices_on_upper_node(3);
407  DenseVector<Number> lower_contact_force_vec(3);
408  DenseVector<Number> upper_contact_force_vec(3);
409 
410  for (unsigned int var=0; var<3; var++)
411  {
412  dof_indices_on_lower_node[var] = lower_node.dof_number(_sys.number(), var, 0);
413  lower_contact_force_vec(var) = -lambda_plus_penalty * contact_force_direction(var);
414 
415  dof_indices_on_upper_node[var] = upper_node.dof_number(_sys.number(), var, 0);
416  upper_contact_force_vec(var) = lambda_plus_penalty * contact_force_direction(var);
417  }
418 
419  if (lambda_plus_penalty > 0.)
420  {
421  if (residual && (_sys.comm().rank() == 0))
422  {
423  residual->add_vector (lower_contact_force_vec, dof_indices_on_lower_node);
424  residual->add_vector (upper_contact_force_vec, dof_indices_on_upper_node);
425  }
426 
427  // Add the Jacobian terms due to the contact forces. The lambda contribution
428  // is not relevant here because it doesn't depend on the solution.
429  if (jacobian && (_sys.comm().rank() == 0))
430  for (unsigned int var=0; var<3; var++)
431  for (unsigned int j=0; j<3; j++)
432  {
433  jacobian->add(dof_indices_on_lower_node[var],
434  dof_indices_on_upper_node[j],
435  _contact_penalty * contact_force_direction(j) * contact_force_direction(var));
436 
437  jacobian->add(dof_indices_on_lower_node[var],
438  dof_indices_on_lower_node[j],
439  -_contact_penalty * contact_force_direction(j) * contact_force_direction(var));
440 
441  jacobian->add(dof_indices_on_upper_node[var],
442  dof_indices_on_lower_node[j],
443  _contact_penalty * contact_force_direction(j) * contact_force_direction(var));
444 
445  jacobian->add(dof_indices_on_upper_node[var],
446  dof_indices_on_upper_node[j],
447  -_contact_penalty * contact_force_direction(j) * contact_force_direction(var));
448  }
449  }
450  else
451  {
452  // We add zeros to the matrix even when lambda_plus_penalty = 0.
453  // We do this because some linear algebra libraries (e.g. PETSc)
454  // will condense out any unused entries from the sparsity pattern,
455  // so adding these zeros in ensures that these entries are not
456  // condensed out.
457  if (jacobian && (_sys.comm().rank() == 0))
458  for (unsigned int var=0; var<3; var++)
459  for (unsigned int j=0; j<3; j++)
460  {
461  jacobian->add(dof_indices_on_lower_node[var],
462  dof_indices_on_upper_node[j],
463  0.);
464 
465  jacobian->add(dof_indices_on_lower_node[var],
466  dof_indices_on_lower_node[j],
467  0.);
468 
469  jacobian->add(dof_indices_on_upper_node[var],
470  dof_indices_on_lower_node[j],
471  0.);
472 
473  jacobian->add(dof_indices_on_upper_node[var],
474  dof_indices_on_upper_node[j],
475  0.);
476  }
477  }
478  }
479 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
A Node is like a Point, but with more information.
Definition: node.h:52
std::map< dof_id_type, Real > _lambda_plus_penalty_values
Store the intermediate values of lambda plus penalty.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
MeshBase & mesh
Real _contact_penalty
Penalize overlapping elements.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
const T & get(const std::string &) const
Definition: parameters.h:430
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
void move_mesh(MeshBase &input_mesh, const NumericVector< Number > &input_solution)
Move the mesh nodes of input_mesh based on the displacement field in input_solution.
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value)=0
Add value to the element (i,j).
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
NonlinearImplicitSystem & _sys
Keep a reference to the NonlinearImplicitSystem.
virtual void zero()=0
Set all entries to 0.
std::map< dof_id_type, Real > _lambdas
Augmented Lagrangian values at each contact node.
Order default_quadrature_order() const
Definition: fe_type.h:332
Defines a dense submatrix for use in Finite Element-type computations.
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:1806
Real elasticity_tensor(Real young_modulus, Real poisson_ratio, unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
unsigned int number() const
Definition: system.h:2006
std::map< dof_id_type, dof_id_type > _contact_node_map
This provides a map between contact nodes.
const Parallel::Communicator & comm() const
This class implements specific orders of Gauss quadrature.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
unsigned int rank() const
Definition: parallel.h:724
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
uint8_t dof_id_type
Definition: id_types.h:64
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
void LinearElasticityWithContact::set_contact_penalty ( Real  contact_penalty_in)

Update the penalty parameter.

Definition at line 40 of file linear_elasticity_with_contact.C.

References _contact_penalty.

41 {
42  _contact_penalty = contact_penalty_in;
43 }
Real _contact_penalty
Penalize overlapping elements.
std::pair< Real, Real > LinearElasticityWithContact::update_lambdas ( )

Update the lambda parameters in the augmented Lagrangian method.

Returns
the largest change in the lambdas, and the largest lambda value.

Definition at line 607 of file linear_elasticity_with_contact.C.

References _lambda_plus_penalty_values, _lambdas, std::abs(), and libMesh::Real.

608 {
609  Real max_delta_lambda = 0.;
610  Real max_new_lambda = 0.;
611 
612  std::map<dof_id_type, Real>::iterator it = _lambdas.begin();
613  std::map<dof_id_type, Real>::iterator it_end = _lambdas.end();
614  for ( ; it != it_end; ++it)
615  {
616  dof_id_type upper_node_id = it->first;
617 
618  std::map<dof_id_type, Real>::iterator new_lambda_it = _lambda_plus_penalty_values.find(upper_node_id);
619  if (new_lambda_it == _lambda_plus_penalty_values.end())
620  libmesh_error_msg("New lambda value not found");
621 
622  Real new_lambda = new_lambda_it->second;
623  Real old_lambda = it->second;
624 
625  it->second = new_lambda;
626 
627  Real delta_lambda = std::abs(new_lambda-old_lambda);
628  if (delta_lambda > max_delta_lambda)
629  max_delta_lambda = delta_lambda;
630 
631  if (std::abs(new_lambda) > max_new_lambda)
632  max_new_lambda = std::abs(new_lambda);
633  }
634 
635  return std::make_pair(max_delta_lambda, max_new_lambda);
636 }
double abs(double a)
std::map< dof_id_type, Real > _lambda_plus_penalty_values
Store the intermediate values of lambda plus penalty.
std::map< dof_id_type, Real > _lambdas
Augmented Lagrangian values at each contact node.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
uint8_t dof_id_type
Definition: id_types.h:64

Member Data Documentation

std::map<dof_id_type, dof_id_type> LinearElasticityWithContact::_contact_node_map
private

This provides a map between contact nodes.

Definition at line 61 of file linear_elasticity_with_contact.h.

Referenced by add_contact_edge_elements(), get_least_and_max_gap_function(), initialize_contact_load_paths(), and residual_and_jacobian().

Real LinearElasticityWithContact::_contact_penalty
private

Penalize overlapping elements.

Definition at line 44 of file linear_elasticity_with_contact.h.

Referenced by get_contact_penalty(), residual_and_jacobian(), and set_contact_penalty().

std::map<dof_id_type, Real> LinearElasticityWithContact::_lambda_plus_penalty_values
private

Store the intermediate values of lambda plus penalty.

The dof IDs refer to the nodes on the upper contact surface.

Definition at line 50 of file linear_elasticity_with_contact.h.

Referenced by residual_and_jacobian(), and update_lambdas().

std::map<dof_id_type, Real> LinearElasticityWithContact::_lambdas
private

Augmented Lagrangian values at each contact node.

The dof IDs refer to the nodes on the upper contact surface.

Definition at line 56 of file linear_elasticity_with_contact.h.

Referenced by initialize_contact_load_paths(), residual_and_jacobian(), and update_lambdas().

NonlinearImplicitSystem& LinearElasticityWithContact::_sys
private

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