libMesh
systems_of_equations_ex7.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 // <h1> Systems Example 7 - Large deformation elasticity (St. Venant-Kirchoff material) </h1>
20 // \author Lorenzo Zanon
21 // \author David Knezevic
22 // \date 2014
23 //
24 // In this example, we consider an elastic cantilever beam modeled as a St. Venant-Kirchoff
25 // material (which is an extension of the linear elastic material model to the nonlinear
26 // regime). The implementation presented here uses NonlinearImplicitSystem.
27 //
28 // We formulate the PDE on the reference geometry (\Omega) as opposed to the deformed
29 // geometry (\Omega^deformed). As a result (e.g. see Ciarlet's 3D elasticity book,
30 // Theorem 2.6-2) the PDE is given as follows:
31 //
32 // \int_\Omega F_im Sigma_mj v_i,j = \int_\Omega f_i v_i + \int_\Gamma g_i v_i ds
33 //
34 // where:
35 // * F is the deformation gradient, F = I + du/dx (x here refers to reference coordinates).
36 // * Sigma is the second Piola-Kirchoff stress, which for the St. Venant Kirchoff model is
37 // given by Sigma_ij = C_ijkl E_kl, where E_kl is the strain,
38 // E_kl = 0.5 * (u_k,l + u_l,k + u_m,k u_m,l).
39 // * f is a body load.
40 // * g is a surface traction on the surface \Gamma.
41 //
42 // In this example we only consider a body load (e.g. gravity), hence we set g = 0.
43 
44 // C++ include files that we need
45 #include <iostream>
46 #include <algorithm>
47 #include <cmath>
48 
49 // Various include files needed for the mesh & solver functionality.
50 #include "libmesh/libmesh.h"
51 #include "libmesh/mesh.h"
52 #include "libmesh/mesh_refinement.h"
53 #include "libmesh/exodusII_io.h"
54 #include "libmesh/equation_systems.h"
55 #include "libmesh/fe.h"
56 #include "libmesh/quadrature_gauss.h"
57 #include "libmesh/dof_map.h"
58 #include "libmesh/sparse_matrix.h"
59 #include "libmesh/numeric_vector.h"
60 #include "libmesh/dense_matrix.h"
61 #include "libmesh/dense_vector.h"
62 #include "libmesh/elem.h"
63 #include "libmesh/string_to_enum.h"
64 #include "libmesh/getpot.h"
65 #include "libmesh/mesh_generation.h"
66 #include "libmesh/dirichlet_boundaries.h"
67 #include "libmesh/zero_function.h"
68 
69 // The nonlinear solver and system we will be using
70 #include "libmesh/nonlinear_solver.h"
71 #include "libmesh/nonlinear_implicit_system.h"
72 
73 #define BOUNDARY_ID_MIN_Z 0
74 #define BOUNDARY_ID_MIN_Y 1
75 #define BOUNDARY_ID_MAX_X 2
76 #define BOUNDARY_ID_MAX_Y 3
77 #define BOUNDARY_ID_MIN_X 4
78 #define BOUNDARY_ID_MAX_Z 5
79 
80 using namespace libMesh;
81 
82 
83 
86 {
87 private:
89 
90 public:
91 
93  es(es_in)
94  {}
95 
99  Real kronecker_delta(unsigned int i,
100  unsigned int j)
101  {
102  return i == j ? 1. : 0.;
103  }
104 
108  Real elasticity_tensor(Real young_modulus,
109  Real poisson_ratio,
110  unsigned int i,
111  unsigned int j,
112  unsigned int k,
113  unsigned int l)
114  {
115  // Define the Lame constants
116  const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
117  const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
118 
119  return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l) +
120  lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l) + kronecker_delta(i,l) * kronecker_delta(j,k));
121  }
122 
123 
127  virtual void jacobian (const NumericVector<Number> & soln,
128  SparseMatrix<Number> & jacobian,
129  NonlinearImplicitSystem & /*sys*/)
130  {
131  const Real young_modulus = es.parameters.get<Real>("young_modulus");
132  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
133 
134  const MeshBase & mesh = es.get_mesh();
135  const unsigned int dim = mesh.mesh_dimension();
136 
137  NonlinearImplicitSystem & system =
138  es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");
139 
140  const unsigned int u_var = system.variable_number ("u");
141 
142  const DofMap & dof_map = system.get_dof_map();
143 
144  FEType fe_type = dof_map.variable_type(u_var);
145  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
146  QGauss qrule (dim, fe_type.default_quadrature_order());
147  fe->attach_quadrature_rule (&qrule);
148 
149  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
150  QGauss qface (dim-1, fe_type.default_quadrature_order());
151  fe_face->attach_quadrature_rule (&qface);
152 
153  const std::vector<Real> & JxW = fe->get_JxW();
154  const std::vector<std::vector<Real>> & phi = fe->get_phi();
155  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
156 
158  DenseSubMatrix<Number> Ke_var[3][3] =
159  {
163  };
164 
165  std::vector<dof_id_type> dof_indices;
166  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
167 
168  jacobian.zero();
169 
170  for (const auto & elem : mesh.active_local_element_ptr_range())
171  {
172  dof_map.dof_indices (elem, dof_indices);
173  for (unsigned int var=0; var<3; var++)
174  dof_map.dof_indices (elem, dof_indices_var[var], var);
175 
176  const unsigned int n_dofs = dof_indices.size();
177  const unsigned int n_var_dofs = dof_indices_var[0].size();
178 
179  fe->reinit (elem);
180 
181  Ke.resize (n_dofs, n_dofs);
182  for (unsigned int var_i=0; var_i<3; var_i++)
183  for (unsigned int var_j=0; var_j<3; var_j++)
184  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
185 
186  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
187  {
188  DenseVector<Number> u_vec(3);
189  DenseMatrix<Number> grad_u(3, 3);
190  for (unsigned int var_i=0; var_i<3; var_i++)
191  {
192  for (unsigned int j=0; j<n_var_dofs; j++)
193  u_vec(var_i) += phi[j][qp]*soln(dof_indices_var[var_i][j]);
194 
195  // Row is variable u1, u2, or u3, column is x, y, or z
196  for (unsigned int var_j=0; var_j<3; var_j++)
197  for (unsigned int j=0; j<n_var_dofs; j++)
198  grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
199  }
200 
201  DenseMatrix<Number> strain_tensor(3, 3);
202  for (unsigned int i=0; i<3; i++)
203  for (unsigned int j=0; j<3; j++)
204  {
205  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));
206 
207  for (unsigned int k=0; k<3; k++)
208  strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
209  }
210 
211  // Define the deformation gradient
212  DenseMatrix<Number> F(3, 3);
213  F = grad_u;
214  for (unsigned int var=0; var<3; var++)
215  F(var, var) += 1.;
216 
217  DenseMatrix<Number> stress_tensor(3, 3);
218 
219  for (unsigned int i=0; i<3; i++)
220  for (unsigned int j=0; j<3; j++)
221  for (unsigned int k=0; k<3; k++)
222  for (unsigned int l=0; l<3; l++)
223  stress_tensor(i,j) +=
224  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k, l);
225 
226  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
227  for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
228  {
229  for (unsigned int i=0; i<3; i++)
230  for (unsigned int j=0; j<3; j++)
231  for (unsigned int m=0; m<3; m++)
232  Ke_var[i][i](dof_i,dof_j) += JxW[qp] *
233  (-dphi[dof_j][qp](m) * stress_tensor(m,j) * dphi[dof_i][qp](j));
234 
235  for (unsigned int i=0; i<3; i++)
236  for (unsigned int j=0; j<3; j++)
237  for (unsigned int k=0; k<3; k++)
238  for (unsigned int l=0; l<3; l++)
239  {
240  Number FxC_ijkl = 0.;
241  for (unsigned int m=0; m<3; m++)
242  FxC_ijkl += F(i,m) * elasticity_tensor(young_modulus, poisson_ratio, m, j, k, l);
243 
244  Ke_var[i][k](dof_i,dof_j) += JxW[qp] *
245  (-0.5 * FxC_ijkl * dphi[dof_j][qp](l) * dphi[dof_i][qp](j));
246 
247  Ke_var[i][l](dof_i,dof_j) += JxW[qp] *
248  (-0.5 * FxC_ijkl * dphi[dof_j][qp](k) * dphi[dof_i][qp](j));
249 
250  for (unsigned int n=0; n<3; n++)
251  Ke_var[i][n](dof_i,dof_j) += JxW[qp] *
252  (-0.5 * FxC_ijkl * (dphi[dof_j][qp](k) * grad_u(n,l) + dphi[dof_j][qp](l) * grad_u(n,k)) * dphi[dof_i][qp](j));
253  }
254  }
255  }
256 
257  dof_map.constrain_element_matrix (Ke, dof_indices);
258  jacobian.add_matrix (Ke, dof_indices);
259  }
260  }
261 
265  virtual void residual (const NumericVector<Number> & soln,
266  NumericVector<Number> & residual,
267  NonlinearImplicitSystem & /*sys*/)
268  {
269  const Real young_modulus = es.parameters.get<Real>("young_modulus");
270  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
271  const Real forcing_magnitude = es.parameters.get<Real>("forcing_magnitude");
272 
273  const MeshBase & mesh = es.get_mesh();
274  const unsigned int dim = mesh.mesh_dimension();
275 
276  NonlinearImplicitSystem & system =
277  es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");
278 
279  const unsigned int u_var = system.variable_number ("u");
280 
281  const DofMap & dof_map = system.get_dof_map();
282 
283  FEType fe_type = dof_map.variable_type(u_var);
284  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
285  QGauss qrule (dim, fe_type.default_quadrature_order());
286  fe->attach_quadrature_rule (&qrule);
287 
288  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
289  QGauss qface (dim-1, fe_type.default_quadrature_order());
290  fe_face->attach_quadrature_rule (&qface);
291 
292  const std::vector<Real> & JxW = fe->get_JxW();
293  const std::vector<std::vector<Real>> & phi = fe->get_phi();
294  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
295 
297 
298  DenseSubVector<Number> Re_var[3] =
302 
303  std::vector<dof_id_type> dof_indices;
304  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
305 
306  residual.zero();
307 
308  for (const auto & elem : mesh.active_local_element_ptr_range())
309  {
310  dof_map.dof_indices (elem, dof_indices);
311  for (unsigned int var=0; var<3; var++)
312  dof_map.dof_indices (elem, dof_indices_var[var], var);
313 
314  const unsigned int n_dofs = dof_indices.size();
315  const unsigned int n_var_dofs = dof_indices_var[0].size();
316 
317  fe->reinit (elem);
318 
319  Re.resize (n_dofs);
320  for (unsigned int var=0; var<3; var++)
321  Re_var[var].reposition (var*n_var_dofs, n_var_dofs);
322 
323  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
324  {
325  DenseVector<Number> u_vec(3);
326  DenseMatrix<Number> grad_u(3, 3);
327  for (unsigned int var_i=0; var_i<3; var_i++)
328  {
329  for (unsigned int j=0; j<n_var_dofs; j++)
330  u_vec(var_i) += phi[j][qp]*soln(dof_indices_var[var_i][j]);
331 
332  // Row is variable u, v, or w column is x, y, or z
333  for (unsigned int var_j=0; var_j<3; var_j++)
334  for (unsigned int j=0; j<n_var_dofs; j++)
335  grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
336  }
337 
338  DenseMatrix<Number> strain_tensor(3, 3);
339  for (unsigned int i=0; i<3; i++)
340  for (unsigned int j=0; j<3; j++)
341  {
342  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));
343 
344  for (unsigned int k=0; k<3; k++)
345  strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
346  }
347 
348  // Define the deformation gradient
349  DenseMatrix<Number> F(3, 3);
350  F = grad_u;
351  for (unsigned int var=0; var<3; var++)
352  F(var, var) += 1.;
353 
354  DenseMatrix<Number> stress_tensor(3, 3);
355 
356  for (unsigned int i=0; i<3; i++)
357  for (unsigned int j=0; j<3; j++)
358  for (unsigned int k=0; k<3; k++)
359  for (unsigned int l=0; l<3; l++)
360  stress_tensor(i,j) +=
361  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k,l);
362 
363  DenseVector<Number> f_vec(3);
364  f_vec(0) = 0.;
365  f_vec(1) = 0.;
366  f_vec(2) = -forcing_magnitude;
367 
368  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
369  for (unsigned int i=0; i<3; i++)
370  {
371  for (unsigned int j=0; j<3; j++)
372  {
373  Number FxStress_ij = 0.;
374  for (unsigned int m=0; m<3; m++)
375  FxStress_ij += F(i,m) * stress_tensor(m,j);
376 
377  Re_var[i](dof_i) += JxW[qp] * (-FxStress_ij * dphi[dof_i][qp](j));
378  }
379 
380  Re_var[i](dof_i) += JxW[qp] * (f_vec(i) * phi[dof_i][qp]);
381  }
382  }
383 
384  dof_map.constrain_element_vector (Re, dof_indices);
385  residual.add_vector (Re, dof_indices);
386  }
387  }
388 
393  {
394  const Real young_modulus = es.parameters.get<Real>("young_modulus");
395  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
396 
397  const MeshBase & mesh = es.get_mesh();
398  const unsigned int dim = mesh.mesh_dimension();
399 
400  NonlinearImplicitSystem & system =
401  es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");
402 
403  unsigned int displacement_vars[3];
404  displacement_vars[0] = system.variable_number ("u");
405  displacement_vars[1] = system.variable_number ("v");
406  displacement_vars[2] = system.variable_number ("w");
407  const unsigned int u_var = system.variable_number ("u");
408 
409  const DofMap & dof_map = system.get_dof_map();
410  FEType fe_type = dof_map.variable_type(u_var);
411  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
412  QGauss qrule (dim, fe_type.default_quadrature_order());
413  fe->attach_quadrature_rule (&qrule);
414 
415  const std::vector<Real> & JxW = fe->get_JxW();
416  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
417 
418  // Also, get a reference to the ExplicitSystem
419  ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
420  const DofMap & stress_dof_map = stress_system.get_dof_map();
421  unsigned int sigma_vars[6];
422  sigma_vars[0] = stress_system.variable_number ("sigma_00");
423  sigma_vars[1] = stress_system.variable_number ("sigma_01");
424  sigma_vars[2] = stress_system.variable_number ("sigma_02");
425  sigma_vars[3] = stress_system.variable_number ("sigma_11");
426  sigma_vars[4] = stress_system.variable_number ("sigma_12");
427  sigma_vars[5] = stress_system.variable_number ("sigma_22");
428 
429  // Storage for the stress dof indices on each element
430  std::vector<std::vector<dof_id_type>> dof_indices_var(system.n_vars());
431  std::vector<dof_id_type> stress_dof_indices_var;
432 
433  // To store the stress tensor on each element
434  DenseMatrix<Number> elem_avg_stress_tensor(3, 3);
435 
436  for (const auto & elem : mesh.active_local_element_ptr_range())
437  {
438  for (unsigned int var=0; var<3; var++)
439  dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
440 
441  const unsigned int n_var_dofs = dof_indices_var[0].size();
442 
443  fe->reinit (elem);
444 
445  // clear the stress tensor
446  elem_avg_stress_tensor.resize(3, 3);
447 
448  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
449  {
450  DenseMatrix<Number> grad_u(3, 3);
451  // Row is variable u1, u2, or u3, column is x, y, or z
452  for (unsigned int var_i=0; var_i<3; var_i++)
453  for (unsigned int var_j=0; var_j<3; var_j++)
454  for (unsigned int j=0; j<n_var_dofs; j++)
455  grad_u(var_i,var_j) += dphi[j][qp](var_j) * system.current_solution(dof_indices_var[var_i][j]);
456 
457  DenseMatrix<Number> strain_tensor(3, 3);
458  for (unsigned int i=0; i<3; i++)
459  for (unsigned int j=0; j<3; j++)
460  {
461  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));
462 
463  for (unsigned int k=0; k<3; k++)
464  strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
465  }
466 
467  // Define the deformation gradient
468  DenseMatrix<Number> F(3, 3);
469  F = grad_u;
470  for (unsigned int var=0; var<3; var++)
471  F(var, var) += 1.;
472 
473  DenseMatrix<Number> stress_tensor(3, 3);
474  for (unsigned int i=0; i<3; i++)
475  for (unsigned int j=0; j<3; j++)
476  for (unsigned int k=0; k<3; k++)
477  for (unsigned int l=0; l<3; l++)
478  stress_tensor(i,j) +=
479  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k, l);
480 
481  // stress_tensor now holds the second Piola-Kirchoff stress (PK2) at point qp.
482  // However, in this example we want to compute the Cauchy stress which is given by
483  // 1/det(F) * F * PK2 * F^T, hence we now apply this transformation.
484  stress_tensor.scale(1./F.det());
485  stress_tensor.left_multiply(F);
486  stress_tensor.right_multiply_transpose(F);
487 
488  // We want to plot the average Cauchy stress on each element, hence
489  // we integrate stress_tensor
490  elem_avg_stress_tensor.add(JxW[qp], stress_tensor);
491  }
492 
493  // Get the average stress per element by dividing by volume
494  elem_avg_stress_tensor.scale(1./elem->volume());
495 
496  // load elem_sigma data into stress_system
497  unsigned int stress_var_index = 0;
498  for (unsigned int i=0; i<3; i++)
499  for (unsigned int j=i; j<3; j++)
500  {
501  stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[stress_var_index]);
502 
503  // We are using CONSTANT MONOMIAL basis functions, hence we only need to get
504  // one dof index per variable
505  dof_id_type dof_index = stress_dof_indices_var[0];
506 
507  if ((stress_system.solution->first_local_index() <= dof_index) &&
508  (dof_index < stress_system.solution->last_local_index()))
509  stress_system.solution->set(dof_index, elem_avg_stress_tensor(i,j));
510 
511  stress_var_index++;
512  }
513  }
514 
515  // Should call close and update when we set vector entries directly
516  stress_system.solution->close();
517  stress_system.update();
518  }
519 
520 };
521 
522 
523 int main (int argc, char ** argv)
524 {
525  LibMeshInit init (argc, argv);
526 
527  // This example requires the PETSc nonlinear solvers
528  libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
529 
530  // We use a 3D domain.
531  libmesh_example_requires(LIBMESH_DIM > 2, "--disable-1D-only --disable-2D-only");
532 
533  GetPot infile("systems_of_equations_ex7.in");
534  const Real x_length = infile("x_length", 0.);
535  const Real y_length = infile("y_length", 0.);
536  const Real z_length = infile("z_length", 0.);
537  const Real n_elem_x = infile("n_elem_x", 0);
538  const Real n_elem_y = infile("n_elem_y", 0);
539  const Real n_elem_z = infile("n_elem_z", 0);
540  const std::string approx_order = infile("approx_order", "FIRST");
541  const std::string fe_family = infile("fe_family", "LAGRANGE");
542 
543  const Real young_modulus = infile("Young_modulus", 1.0);
544  const Real poisson_ratio = infile("poisson_ratio", 0.3);
545  const Real forcing_magnitude = infile("forcing_magnitude", 0.001);
546 
547  const Real nonlinear_abs_tol = infile("nonlinear_abs_tol", 1.e-8);
548  const Real nonlinear_rel_tol = infile("nonlinear_rel_tol", 1.e-8);
549  const unsigned int nonlinear_max_its = infile("nonlinear_max_its", 50);
550 
551  const unsigned int n_solves = infile("n_solves", 10);
552  const Real force_scaling = infile("force_scaling", 5.0);
553 
554  Mesh mesh(init.comm());
555 
557  n_elem_x,
558  n_elem_y,
559  n_elem_z,
560  0., x_length,
561  0., y_length,
562  0., z_length,
563  HEX27);
564 
565  mesh.print_info();
566 
567  EquationSystems equation_systems (mesh);
568  LargeDeformationElasticity lde(equation_systems);
569 
570  NonlinearImplicitSystem & system =
571  equation_systems.add_system<NonlinearImplicitSystem> ("NonlinearElasticity");
572 
573  unsigned int u_var =
574  system.add_variable("u",
575  Utility::string_to_enum<Order> (approx_order),
576  Utility::string_to_enum<FEFamily>(fe_family));
577 
578  unsigned int v_var =
579  system.add_variable("v",
580  Utility::string_to_enum<Order> (approx_order),
581  Utility::string_to_enum<FEFamily>(fe_family));
582 
583  unsigned int w_var =
584  system.add_variable("w",
585  Utility::string_to_enum<Order> (approx_order),
586  Utility::string_to_enum<FEFamily>(fe_family));
587 
588  // Also, initialize an ExplicitSystem to store stresses
589  ExplicitSystem & stress_system =
590  equation_systems.add_system<ExplicitSystem> ("StressSystem");
591  stress_system.add_variable("sigma_00", CONSTANT, MONOMIAL);
592  stress_system.add_variable("sigma_01", CONSTANT, MONOMIAL);
593  stress_system.add_variable("sigma_02", CONSTANT, MONOMIAL);
594  stress_system.add_variable("sigma_11", CONSTANT, MONOMIAL);
595  stress_system.add_variable("sigma_12", CONSTANT, MONOMIAL);
596  stress_system.add_variable("sigma_22", CONSTANT, MONOMIAL);
597 
598  equation_systems.parameters.set<Real> ("nonlinear solver absolute residual tolerance") = nonlinear_abs_tol;
599  equation_systems.parameters.set<Real> ("nonlinear solver relative residual tolerance") = nonlinear_rel_tol;
600  equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = nonlinear_max_its;
601 
602  system.nonlinear_solver->residual_object = &lde;
603  system.nonlinear_solver->jacobian_object = &lde;
604 
605  equation_systems.parameters.set<Real>("young_modulus") = young_modulus;
606  equation_systems.parameters.set<Real>("poisson_ratio") = poisson_ratio;
607  equation_systems.parameters.set<Real>("forcing_magnitude") = forcing_magnitude;
608 
609  // Attach Dirichlet boundary conditions
610  std::set<boundary_id_type> clamped_boundaries;
611  clamped_boundaries.insert(BOUNDARY_ID_MIN_X);
612 
613  std::vector<unsigned int> uvw;
614  uvw.push_back(u_var);
615  uvw.push_back(v_var);
616  uvw.push_back(w_var);
617 
619 
620  // Most DirichletBoundary users will want to supply a "locally
621  // indexed" functor
623  (DirichletBoundary (clamped_boundaries, uvw, zero,
625 
626  equation_systems.init();
627  equation_systems.print_info();
628 
629  // Provide a loop here so that we can do a sequence of solves
630  // where solve n gives a good starting guess for solve n+1.
631  // This "continuation" approach is helpful for solving for
632  // large values of "forcing_magnitude".
633  // Set n_solves and force_scaling in nonlinear_elasticity.in.
634  for (unsigned int count=0; count<n_solves; count++)
635  {
636  Real previous_forcing_magnitude = equation_systems.parameters.get<Real>("forcing_magnitude");
637  equation_systems.parameters.set<Real>("forcing_magnitude") = previous_forcing_magnitude*force_scaling;
638 
639  libMesh::out << "Performing solve "
640  << count
641  << ", forcing_magnitude: "
642  << equation_systems.parameters.get<Real>("forcing_magnitude")
643  << std::endl;
644 
645  system.solve();
646 
647  libMesh::out << "System solved at nonlinear iteration "
648  << system.n_nonlinear_iterations()
649  << " , final nonlinear residual norm: "
650  << system.final_nonlinear_residual()
651  << std::endl
652  << std::endl;
653 
654  libMesh::out << "Computing stresses..." << std::endl;
655 
656  lde.compute_stresses();
657 
658 #ifdef LIBMESH_HAVE_EXODUS_API
659  std::stringstream filename;
660  filename << "solution_" << count << ".exo";
661  ExodusII_IO (mesh).write_equation_systems(filename.str(), equation_systems);
662 #endif
663  }
664 
665  return 0;
666 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This is the EquationSystems class.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
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]...
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
void constrain_element_vector(DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
Constrains the element vector.
Definition: dof_map.h:1802
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Definition: assembly.C:38
MeshBase & mesh
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
const T & get(const std::string &) const
Definition: parameters.h:430
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
Real kronecker_delta(unsigned int i, unsigned int j)
Kronecker delta function.
const Number zero
.
Definition: libmesh.h:178
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
virtual void residual(const NumericVector< Number > &soln, NumericVector< Number > &residual, NonlinearImplicitSystem &)
Evaluate the residual of the nonlinear system.
SolverPackage default_solver_package()
Definition: libmesh.C:995
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
Abstract base class to be used to calculate the residual of a nonlinear system.
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.
UniquePtr< NonlinearSolver< Number > > nonlinear_solver
The NonlinearSolver defines the default interface used to solve the nonlinear_implicit system...
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual void zero()=0
Set all entries to 0.
Abstract base class to be used to calculate the Jacobian of a nonlinear system.
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
Order default_quadrature_order() const
Definition: fe_type.h:332
int main(int argc, char **argv)
boostcopy::enable_if_c< ScalarTraits< T2 >::value, void >::type add(const T2 factor, const DenseMatrix< T3 > &mat)
Adds factor times mat to this matrix.
Definition: dense_matrix.h:883
Defines a dense submatrix for use in Finite Element-type computations.
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
This class provides a specific system class.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
Real kronecker_delta(unsigned int i, unsigned int j)
Definition: assembly.C:32
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void solve() libmesh_override
Assembles & solves the nonlinear system R(x) = 0.
virtual void left_multiply(const DenseMatrixBase< T > &M2) libmesh_override
Performs the operation: (*this) <- M2 * (*this)
Definition: dense_matrix.C:37
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
void right_multiply_transpose(const DenseMatrix< T > &A)
Right multiplies by the transpose of the matrix A.
Definition: dense_matrix.C:257
void scale(const T factor)
Multiplies every element in the matrix by factor.
Definition: dense_matrix.h:851
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
void constrain_element_matrix(DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix.
Definition: dof_map.h:1793
virtual void jacobian(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, NonlinearImplicitSystem &)
Evaluate the Jacobian of the nonlinear system.
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
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.
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
LargeDeformationElasticity(EquationSystems &es_in)
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
virtual void init()
Initialize all the systems.
void compute_stresses()
Compute the Cauchy stress for the current solution.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
virtual unsigned int size() const libmesh_override
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
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
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.