libMesh
solid_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // \author Robert Weidlich
21 // \date Copyright 2012
22 
23 #include "libmesh/boundary_info.h"
24 #include "libmesh/diff_solver.h"
25 #include "libmesh/dof_map.h"
26 #include "libmesh/equation_systems.h"
27 #include "libmesh/fe_base.h"
28 #include "libmesh/fem_context.h"
29 #include "libmesh/getpot.h"
30 #include "libmesh/mesh.h"
31 #include "libmesh/newton_solver.h"
32 #include "libmesh/numeric_vector.h"
33 #include "libmesh/quadrature.h"
34 #include "libmesh/sparse_matrix.h"
35 #include "libmesh/steady_solver.h"
36 #include "libmesh/transient_system.h"
37 #include "libmesh/node.h"
38 #include "libmesh/elem.h"
39 
40 #include "nonlinear_neohooke_cc.h"
41 #include "solid_system.h"
42 
43 // Solaris Studio has no NAN
44 #ifdef __SUNPRO_CC
45 #define NAN (1.0/0.0)
46 #endif
47 
48 // Bring in everything from the libMesh namespace
49 using namespace libMesh;
50 
52  const std::string & name_in,
53  const unsigned int number_in) :
54  FEMSystem(es, name_in, number_in)
55 {
56 
57  // Add a time solver. We are just looking at a steady state problem here.
59 }
60 
62 {
63  System & aux_sys = this->get_equation_systems().get_system("auxiliary");
64 
65  const unsigned int dim = this->get_mesh().mesh_dimension();
66 
67  // Loop over all nodes and copy the location from the current system to
68  // the auxiliary system.
69  for (const auto & node : this->get_mesh().local_node_ptr_range())
70  for (unsigned int d = 0; d < dim; ++d)
71  {
72  unsigned int source_dof = node->dof_number(this->number(), var[d], 0);
73  unsigned int dest_dof = node->dof_number(aux_sys.number(), undefo_var[d], 0);
74  Number value = this->current_local_solution->el(source_dof);
75  aux_sys.current_local_solution->set(dest_dof, value);
76  }
77 }
78 
80 {
81  const unsigned int dim = this->get_mesh().mesh_dimension();
82 
83  // Get the default order of the used elements. Assumption:
84  // Just one type of elements in the mesh.
85  Order order = (*(this->get_mesh().elements_begin()))->default_order();
86 
87  // Add the node positions as primary variables.
88  var[0] = this->add_variable("x", order);
89  var[1] = this->add_variable("y", order);
90  if (dim == 3)
91  var[2] = this->add_variable("z", order);
92  else
93  var[2] = var[1];
94 
95  // Add variables for storing the initial mesh to an auxiliary system.
96  System & aux_sys = this->get_equation_systems().get_system("auxiliary");
97  undefo_var[0] = aux_sys.add_variable("undefo_x", order);
98  undefo_var[1] = aux_sys.add_variable("undefo_y", order);
99  undefo_var[2] = aux_sys.add_variable("undefo_z", order);
100 
101  // Set the time stepping options
102  this->deltat = args("schedule/dt", 0.2);
103 
104  // Do the parent's initialization after variables are defined
105  FEMSystem::init_data();
106 
109  //this->time_evolving(var[0]);
110  //this->time_evolving(var[1]);
111  //if (dim == 3)
112  //this->time_evolving(var[2]);
113 
114  // Tell the system which variables are containing the node positions
115  set_mesh_system(this);
116 
117  this->set_mesh_x_var(var[0]);
118  this->set_mesh_y_var(var[1]);
119  if (dim == 3)
120  this->set_mesh_z_var(var[2]);
121 
122  // Fill the variables with the position of the nodes
123  this->mesh_position_get();
124 
125  System::reinit();
126 
127  // Set some options for the DiffSolver
128  DiffSolver & solver = *(this->time_solver->diff_solver().get());
129  solver.quiet = args("solver/quiet", false);
130  solver.max_nonlinear_iterations = args("solver/nonlinear/max_nonlinear_iterations", 100);
131  solver.relative_step_tolerance = args("solver/nonlinear/relative_step_tolerance", 1.e-3);
132  solver.relative_residual_tolerance = args("solver/nonlinear/relative_residual_tolerance", 1.e-8);
133  solver.absolute_residual_tolerance = args("solver/nonlinear/absolute_residual_tolerance", 1.e-8);
134  solver.verbose = !args("solver/quiet", false);
135 
136  ((NewtonSolver &) solver).require_residual_reduction = args("solver/nonlinear/require_reduction", false);
137 
138  // And the linear solver options
139  solver.max_linear_iterations = args("max_linear_iterations", 50000);
140  solver.initial_linear_tolerance = args("initial_linear_tolerance", 1.e-3);
141 }
142 
144 {
145  System::update();
146  this->mesh_position_set();
147 }
148 
150 {
151  FEMContext & c = cast_ref<FEMContext &>(context);
152 
153  // Pre-request all the data needed
154  FEBase * elem_fe = libmesh_nullptr;
155  c.get_element_fe(0, elem_fe);
156 
157  elem_fe->get_JxW();
158  elem_fe->get_phi();
159  elem_fe->get_dphi();
160  elem_fe->get_xyz();
161 
162  FEBase * side_fe = libmesh_nullptr;
163  c.get_side_fe(0, side_fe);
164 
165  side_fe->get_JxW();
166  side_fe->get_phi();
167  side_fe->get_xyz();
168 }
169 
174 bool SolidSystem::element_time_derivative(bool request_jacobian,
175  DiffContext & context)
176 {
177  FEMContext & c = cast_ref<FEMContext &>(context);
178 
179  // First we get some references to cell-specific data that
180  // will be used to assemble the linear system.
181  FEBase * elem_fe = libmesh_nullptr;
182  c.get_element_fe(0, elem_fe);
183 
184  // Element Jacobian * quadrature weights for interior integration
185  const std::vector<Real> & JxW = elem_fe->get_JxW();
186 
187  // Element basis functions
188  const std::vector<std::vector<RealGradient>> & dphi = elem_fe->get_dphi();
189 
190  // Dimension of the mesh
191  const unsigned int dim = this->get_mesh().mesh_dimension();
192 
193  // The number of local degrees of freedom in each variable
194  const unsigned int n_u_dofs = c.get_dof_indices(var[0]).size();
195  libmesh_assert(n_u_dofs == c.get_dof_indices(var[1]).size());
196  if (dim == 3)
197  libmesh_assert(n_u_dofs == c.get_dof_indices(var[2]).size());
198 
199  unsigned int n_qpoints = c.get_element_qrule().n_points();
200 
201  // Some matrices and vectors for storing the results of the constitutive
202  // law
203  DenseMatrix<Real> stiff;
204  DenseVector<Real> res;
205  VectorValue<Gradient> grad_u;
206 
207  // Instantiate the constitutive law
208  // Just calculate jacobian contribution when we need to
209  NonlinearNeoHookeCurrentConfig material(dphi, args, request_jacobian);
210 
211  // Get a reference to the auxiliary system
213  TransientExplicitSystem>("auxiliary");
214  std::vector<dof_id_type> undefo_index;
215 
216  // Assume symmetry of local stiffness matrices
217  bool use_symmetry = args("assembly/use_symmetry", false);
218 
219  // Now we will build the element Jacobian and residual.
220  // This must be calculated at each quadrature point by
221  // summing the solution degree-of-freedom values by
222  // the appropriate weight functions.
223  // This class just takes care of the assembly. The matrix of
224  // the jacobian and the residual vector are provided by the
225  // constitutive formulation.
226 
227  for (unsigned int qp = 0; qp != n_qpoints; qp++)
228  {
229  // Compute the displacement gradient
230  grad_u(0) = grad_u(1) = grad_u(2) = 0;
231  for (unsigned int d = 0; d < dim; ++d)
232  {
233  std::vector<Number> u_undefo;
234  aux_system.get_dof_map().dof_indices(&c.get_elem(), undefo_index, undefo_var[d]);
235  aux_system.current_local_solution->get(undefo_index, u_undefo);
236  for (unsigned int l = 0; l != n_u_dofs; l++)
237  grad_u(d).add_scaled(dphi[l][qp], u_undefo[l]); // u_current(l)); // -
238  }
239 
240  // initialize the constitutive formulation with the current displacement
241  // gradient
242  material.init_for_qp(grad_u, qp);
243 
244  // Acquire, scale and assemble residual and stiffness
245  for (unsigned int i = 0; i < n_u_dofs; i++)
246  {
247  res.resize(dim);
248  material.get_residual(res, i);
249  res.scale(JxW[qp]);
250  for (unsigned int ii = 0; ii < dim; ++ii)
251  (c.get_elem_residual(ii))(i) += res(ii);
252 
253  if (request_jacobian && c.elem_solution_derivative)
254  {
256  for (unsigned int j = (use_symmetry ? i : 0); j < n_u_dofs; j++)
257  {
258  material.get_linearized_stiffness(stiff, i, j);
259  stiff.scale(JxW[qp]);
260  for (unsigned int ii = 0; ii < dim; ++ii)
261  {
262  for (unsigned int jj = 0; jj < dim; ++jj)
263  {
264  (c.get_elem_jacobian(ii,jj))(i, j) += stiff(ii, jj);
265  if (use_symmetry && i != j)
266  (c.get_elem_jacobian(ii,jj))(j, i) += stiff(jj, ii);
267  }
268  }
269  }
270  }
271  }
272  } // end of the quadrature point qp-loop
273 
274  return request_jacobian;
275 }
276 
277 bool SolidSystem::side_time_derivative(bool request_jacobian,
278  DiffContext & context)
279 {
280  FEMContext & c = cast_ref<FEMContext &>(context);
281 
282  // Apply displacement boundary conditions with penalty method
283 
284  // Get the current load step
285  Real ratio = this->get_equation_systems().parameters.get<Real>("progress") + 0.001;
286 
287  // The BC are stored in the simulation parameters as array containing sequences of
288  // four numbers: Id of the side for the displacements and three values describing the
289  // displacement. E.g.: bc/displacement = '5 nan nan -1.0'. This will move all nodes of
290  // side 5 about 1.0 units down the z-axis while leaving all other directions unrestricted
291 
292  // Get number of BCs to enforce
293  unsigned int num_bc = args.vector_variable_size("bc/displacement");
294  if (num_bc % 4 != 0)
295  libmesh_error_msg("ERROR, Odd number of values in displacement boundary condition.");
296  num_bc /= 4;
297 
298  // Loop over all BCs
299  for (unsigned int nbc = 0; nbc < num_bc; nbc++)
300  {
301  // Get IDs of the side for this BC
302  short int positive_boundary_id = args("bc/displacement", 1, nbc * 4);
303 
304  // The current side may not be on the boundary to be restricted
305  if (!this->get_mesh().get_boundary_info().has_boundary_id
306  (&c.get_elem(),c.get_side(),positive_boundary_id))
307  continue;
308 
309  // Read values from configuration file
310  Point diff_value;
311  for (unsigned int d = 0; d < c.get_dim(); ++d)
312  diff_value(d) = args("bc/displacement", NAN, nbc * 4 + 1 + d);
313 
314  // Scale according to current load step
315  diff_value *= ratio;
316 
317  Real penalty_number = args("bc/displacement_penalty", 1e7);
318 
319  FEBase * fe = libmesh_nullptr;
320  c.get_side_fe(0, fe);
321 
322  const std::vector<std::vector<Real>> & phi = fe->get_phi();
323  const std::vector<Real> & JxW = fe->get_JxW();
324  const std::vector<Point> & coords = fe->get_xyz();
325 
326  unsigned int n_x_dofs = c.get_dof_indices(this->var[0]).size();
327 
328  // get mappings for dofs for auxiliary system for original mesh positions
329  const System & auxsys = this->get_equation_systems().get_system("auxiliary");
330  const DofMap & auxmap = auxsys.get_dof_map();
331  std::vector<dof_id_type> undefo_dofs[3];
332  for (unsigned int d = 0; d < c.get_dim(); ++d)
333  auxmap.dof_indices(&c.get_elem(), undefo_dofs[d], undefo_var[d]);
334 
335  for (unsigned int qp = 0; qp < c.get_side_qrule().n_points(); ++qp)
336  {
337  // calculate coordinates of qp on undeformed mesh
338  Point orig_point;
339  for (unsigned int i = 0; i < n_x_dofs; ++i)
340  {
341  for (unsigned int d = 0; d < c.get_dim(); ++d)
342  {
343  Number orig_val = auxsys.current_solution(undefo_dofs[d][i]);
344 
345 #if LIBMESH_USE_COMPLEX_NUMBERS
346  orig_point(d) += phi[i][qp] * orig_val.real();
347 #else
348  orig_point(d) += phi[i][qp] * orig_val;
349 #endif
350  }
351  }
352 
353  // Calculate displacement to be enforced.
354  Point diff = coords[qp] - orig_point - diff_value;
355 
356  // Assemble
357  for (unsigned int i = 0; i < n_x_dofs; ++i)
358  {
359  for (unsigned int d1 = 0; d1 < c.get_dim(); ++d1)
360  {
361  if (libmesh_isnan(diff(d1)))
362  continue;
363  Real val = JxW[qp] * phi[i][qp] * diff(d1) * penalty_number;
364  (c.get_elem_residual(var[d1]))(i) += val;
365  }
366  if (request_jacobian)
367  {
368  for (unsigned int j = 0; j < n_x_dofs; ++j)
369  {
370  for (unsigned int d1 = 0; d1 < c.get_dim(); ++d1)
371  {
372  if (libmesh_isnan(diff(d1)))
373  continue;
374  Real val = JxW[qp] * phi[i][qp] * phi[j][qp] * penalty_number;
375  (c.get_elem_jacobian(var[d1],var[d1]))(i, j) += val;
376  }
377  }
378  }
379  }
380  }
381  }
382 
383  return request_jacobian;
384 }
unsigned char get_side() const
Accessor for current side of Elem object.
Definition: fem_context.h:885
UniquePtr< TimeSolver > time_solver
A pointer to the solver object we&#39;re going to use.
Definition: diff_system.h:221
This is the EquationSystems class.
void init_for_qp(VectorValue< Gradient > &grad_u, unsigned int qp)
Initialize the class for the given displacement gradient at the specified quadrature point...
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Compute contribution from internal forces in elem_residual and contribution from linearized internal ...
Definition: solid_system.C:174
virtual bool side_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on side of elem to elem_residual.
Definition: solid_system.C:277
This class implements a constitutive formulation for an Neo-Hookean elastic solid in terms of the cur...
bool quiet
The DiffSolver should not print anything to libMesh::out unless quiet is set to false; default is tru...
Definition: diff_solver.h:162
void add_scaled(const TypeVector< T2 > &, const T)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:624
unsigned int dim
void get_residual(DenseVector< Real > &residuum, unsigned int &i)
Return the residual vector for the current state.
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:772
unsigned int max_nonlinear_iterations
The DiffSolver should exit in failure if max_nonlinear_iterations is exceeded and continue_after_max_...
Definition: diff_solver.h:156
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
Real absolute_residual_tolerance
The DiffSolver should exit after the residual is reduced to either less than absolute_residual_tolera...
Definition: diff_solver.h:191
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
const class libmesh_nullptr_t libmesh_nullptr
unsigned int var[3]
Definition: solid_system.h:74
This class provides a specific system class.
const T & get(const std::string &) const
Definition: parameters.h:430
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The libMesh namespace provides an interface to certain functionality in the library.
Real initial_linear_tolerance
Any required linear solves will at first be done with this tolerance; the DiffSolver may tighten the ...
Definition: diff_solver.h:210
void scale(const T factor)
Multiplies every element in the vector by factor.
Definition: dense_vector.h:407
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
virtual void init_context(DiffContext &context)
Definition: solid_system.C:149
This class provides a specific system class.
Definition: fem_system.h:53
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
This class defines a solver which uses the default libMesh linear solver in a quasiNewton method to h...
Definition: newton_solver.h:46
virtual void set_mesh_x_var(unsigned int var)
Tells the DifferentiablePhysics that variable var from the mesh system should be used to update the x...
Definition: diff_physics.h:599
libmesh_assert(j)
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: solid_system.C:79
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
unsigned char get_dim() const
Accessor for cached mesh dimension.
Definition: fem_context.h:899
PetscDiffSolver & solver
virtual element_iterator elements_begin()=0
Iterate over all the elements in the Mesh.
This is a generic class that defines a solver to handle ImplicitSystem classes, including NonlinearIm...
Definition: diff_solver.h:70
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
GetPot args
Definition: solid_system.h:59
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
const MeshBase & get_mesh() const
Definition: system.h:2014
const DofMap & get_dof_map() const
Definition: system.h:2030
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
void get_linearized_stiffness(DenseMatrix< Real > &stiffness, unsigned int &i, unsigned int &j)
Return the stiffness matrix for the current state.
virtual void set_mesh_y_var(unsigned int var)
Tells the DifferentiablePhysics that variable var from the mesh system should be used to update the y...
Definition: diff_physics.h:607
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: solid_system.C:143
Real deltat
For time-dependent problems, this is the amount delta t to advance the solution in time...
Definition: diff_system.h:248
unsigned int max_linear_iterations
Each linear solver step should exit after max_linear_iterations is exceeded.
Definition: diff_solver.h:148
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
Real elem_solution_derivative
The derivative of elem_solution with respect to the current nonlinear solution.
Definition: diff_context.h:483
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
This class implements a TimeSolver which does a single solve of the steady state problem.
Definition: steady_solver.h:47
SolidSystem(EquationSystems &es, const std::string &name, const unsigned int number)
Definition: solid_system.C:51
unsigned int undefo_var[3]
Definition: solid_system.h:78
void scale(const T factor)
Multiplies every element in the matrix by factor.
Definition: dense_matrix.h:851
bool verbose
The DiffSolver may print a lot more to libMesh::out if verbose is set to true; default is false...
Definition: diff_solver.h:168
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
unsigned int number() const
Definition: system.h:2006
bool libmesh_isnan(float a)
virtual void set_mesh_z_var(unsigned int var)
Tells the DifferentiablePhysics that variable var from the mesh system should be used to update the z...
Definition: diff_physics.h:615
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void mesh_position_set()
Tells the FEMSystem to set the mesh nodal coordinates which should correspond to degree of freedom co...
Definition: fem_system.C:1067
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
static const bool value
Definition: xdr_io.C:108
void save_initial_mesh()
Definition: solid_system.C:61
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
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
Defines a dense vector for use in Finite Element-type computations.
void mesh_position_get()
Tells the FEMSystem to set the degree of freedom coefficients which should correspond to mesh nodal c...
Definition: fem_system.C:1392
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:299
virtual void set_mesh_system(System *sys)
Tells the DifferentiablePhysics that system sys contains the isoparametric Lagrangian variables which...
Definition: diff_physics.h:579
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int n_points() const
Definition: quadrature.h:113
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
This class forms the foundation from which generic finite elements may be derived.
Real relative_residual_tolerance
Definition: diff_solver.h:192
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