libMesh
systems_of_equations_ex6.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 // <h1> Systems Example 6 - 3D Linear Elastic Cantilever </h1>
21 // \author David Knezevic
22 // \date 2012
23 //
24 // This is a 3D version of systems_of_equations_ex4. The weak form PDE for
25 // equilibrium elasticity is:
26 //
27 // \int_\Omega Sigma_ij v_i,j = \int_\Omega f_i v_i + \int_\Gamma g_i v_i ds,
28 //
29 // for all admissible test functions v, where:
30 // * Sigma is the stress tensor, which for linear elasticity is
31 // given by Sigma_ij = C_ijkl u_k,l.
32 // * f is a body load.
33 // * g is a surface traction on the surface \Gamma.
34 
35 
36 // C++ include files that we need
37 #include <iostream>
38 #include <algorithm>
39 #include <math.h>
40 
41 // libMesh includes
42 #include "libmesh/libmesh_config.h"
43 #include "libmesh/libmesh.h"
44 #include "libmesh/mesh.h"
45 #include "libmesh/mesh_generation.h"
46 #include "libmesh/exodusII_io.h"
47 #include "libmesh/gnuplot_io.h"
48 #include "libmesh/linear_implicit_system.h"
49 #include "libmesh/equation_systems.h"
50 #include "libmesh/fe.h"
51 #include "libmesh/quadrature_gauss.h"
52 #include "libmesh/dof_map.h"
53 #include "libmesh/sparse_matrix.h"
54 #include "libmesh/numeric_vector.h"
55 #include "libmesh/dense_matrix.h"
56 #include "libmesh/dense_submatrix.h"
57 #include "libmesh/dense_vector.h"
58 #include "libmesh/dense_subvector.h"
59 #include "libmesh/perf_log.h"
60 #include "libmesh/elem.h"
61 #include "libmesh/boundary_info.h"
62 #include "libmesh/zero_function.h"
63 #include "libmesh/dirichlet_boundaries.h"
64 #include "libmesh/string_to_enum.h"
65 #include "libmesh/getpot.h"
66 #include "libmesh/solver_configuration.h"
67 #include "libmesh/petsc_linear_solver.h"
68 #include "libmesh/petsc_macro.h"
69 
70 #define x_scaling 1.3
71 
72 // boundary IDs
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 #define NODE_BOUNDARY_ID 10
80 #define EDGE_BOUNDARY_ID 20
81 
82 // Bring in everything from the libMesh namespace
83 using namespace libMesh;
84 
85 #ifdef LIBMESH_HAVE_PETSC
86 // This class allows us to set the solver and preconditioner
87 // to be appropriate for linear elasticity.
89 {
90 public:
91 
93  _petsc_linear_solver(petsc_linear_solver)
94  {
95  }
96 
97  virtual void configure_solver()
98  {
99  PetscErrorCode ierr = 0;
100  ierr = KSPSetType (_petsc_linear_solver.ksp(), const_cast<KSPType>(KSPCG));
101  libmesh_assert(ierr == 0);
102 
103  ierr = PCSetType (_petsc_linear_solver.pc(), const_cast<PCType>(PCBJACOBI));
104  libmesh_assert(ierr == 0);
105  }
106 
107  // The linear solver object that we are configuring
109 
110 };
111 #endif
112 
114 {
115 private:
117 
118 public:
119 
121  es(es_in)
122  {}
123 
127  Real kronecker_delta(unsigned int i,
128  unsigned int j)
129  {
130  return i == j ? 1. : 0.;
131  }
132 
136  Real elasticity_tensor(unsigned int i,
137  unsigned int j,
138  unsigned int k,
139  unsigned int l)
140  {
141  // Hard code material parameters for the sake of simplicity
142  const Real poisson_ratio = 0.3;
143  const Real young_modulus = 1.;
144 
145  // Define the Lame constants
146  const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
147  const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
148 
149  return lambda_1 * kronecker_delta(i, j) * kronecker_delta(k, l) +
150  lambda_2 * (kronecker_delta(i, k) * kronecker_delta(j, l) + kronecker_delta(i, l) * kronecker_delta(j, k));
151  }
152 
156  void assemble()
157  {
158  const MeshBase & mesh = es.get_mesh();
159 
160  const unsigned int dim = mesh.mesh_dimension();
161 
162  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Elasticity");
163 
164  const unsigned int u_var = system.variable_number ("u");
165 
166  const DofMap & dof_map = system.get_dof_map();
167  FEType fe_type = dof_map.variable_type(u_var);
168  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
169  QGauss qrule (dim, fe_type.default_quadrature_order());
170  fe->attach_quadrature_rule (&qrule);
171 
172  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
173  QGauss qface(dim-1, fe_type.default_quadrature_order());
174  fe_face->attach_quadrature_rule (&qface);
175 
176  const std::vector<Real> & JxW = fe->get_JxW();
177  const std::vector<std::vector<Real>> & phi = fe->get_phi();
178  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
179 
181  DenseSubMatrix<Number> Ke_var[3][3] =
182  {
186  };
187 
189 
190  DenseSubVector<Number> Fe_var[3] =
194 
195  std::vector<dof_id_type> dof_indices;
196  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
197 
198  for (const auto & elem : mesh.active_local_element_ptr_range())
199  {
200  dof_map.dof_indices (elem, dof_indices);
201  for (unsigned int var=0; var<3; var++)
202  dof_map.dof_indices (elem, dof_indices_var[var], var);
203 
204  const unsigned int n_dofs = dof_indices.size();
205  const unsigned int n_var_dofs = dof_indices_var[0].size();
206 
207  fe->reinit (elem);
208 
209  Ke.resize (n_dofs, n_dofs);
210  for (unsigned int var_i=0; var_i<3; var_i++)
211  for (unsigned int var_j=0; var_j<3; var_j++)
212  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
213 
214  Fe.resize (n_dofs);
215  for (unsigned int var=0; var<3; var++)
216  Fe_var[var].reposition (var*n_var_dofs, n_var_dofs);
217 
218  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
219  {
220  // assemble \int_Omega C_ijkl u_k,l v_i,j \dx
221  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
222  for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
223  for (unsigned int i=0; i<3; i++)
224  for (unsigned int j=0; j<3; j++)
225  for (unsigned int k=0; k<3; k++)
226  for (unsigned int l=0; l<3; l++)
227  Ke_var[i][k](dof_i,dof_j) +=
228  JxW[qp] * elasticity_tensor(i,j,k,l) * dphi[dof_j][qp](l) * dphi[dof_i][qp](j);
229 
230  // assemble \int_Omega f_i v_i \dx
231  DenseVector<Number> f_vec(3);
232  f_vec(0) = 0.;
233  f_vec(1) = 0.;
234  f_vec(2) = -1.;
235  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
236  for (unsigned int i=0; i<3; i++)
237  Fe_var[i](dof_i) += JxW[qp] * (f_vec(i) * phi[dof_i][qp]);
238  }
239 
240  // assemble \int_\Gamma g_i v_i \ds
241  DenseVector<Number> g_vec(3);
242  g_vec(0) = 0.;
243  g_vec(1) = 0.;
244  g_vec(2) = -1.;
245  {
246  for (auto side : elem->side_index_range())
247  if (elem->neighbor_ptr(side) == libmesh_nullptr)
248  {
249  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
250  const std::vector<Real> & JxW_face = fe_face->get_JxW();
251 
252  fe_face->reinit(elem, side);
253 
254  // Apply a traction
255  for (unsigned int qp=0; qp<qface.n_points(); qp++)
256  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
257  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
258  for (unsigned int i=0; i<3; i++)
259  Fe_var[i](dof_i) += JxW_face[qp] * (g_vec(i) * phi_face[dof_i][qp]);
260  }
261  }
262 
263  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
264 
265  system.matrix->add_matrix (Ke, dof_indices);
266  system.rhs->add_vector (Fe, dof_indices);
267  }
268  }
269 
270  // Post-process the solution to compute stresses
272  {
273  const MeshBase & mesh = es.get_mesh();
274  const unsigned int dim = mesh.mesh_dimension();
275 
276  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Elasticity");
277 
278  unsigned int displacement_vars[3];
279  displacement_vars[0] = system.variable_number ("u");
280  displacement_vars[1] = system.variable_number ("v");
281  displacement_vars[2] = system.variable_number ("w");
282  const unsigned int u_var = system.variable_number ("u");
283 
284  const DofMap & dof_map = system.get_dof_map();
285  FEType fe_type = dof_map.variable_type(u_var);
286  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
287  QGauss qrule (dim, fe_type.default_quadrature_order());
288  fe->attach_quadrature_rule (&qrule);
289 
290  const std::vector<Real> & JxW = fe->get_JxW();
291  const std::vector<std::vector<Real>> & phi = fe->get_phi();
292  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
293 
294  // Also, get a reference to the ExplicitSystem
295  ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
296  const DofMap & stress_dof_map = stress_system.get_dof_map();
297  unsigned int sigma_vars[6];
298  sigma_vars[0] = stress_system.variable_number ("sigma_00");
299  sigma_vars[1] = stress_system.variable_number ("sigma_01");
300  sigma_vars[2] = stress_system.variable_number ("sigma_02");
301  sigma_vars[3] = stress_system.variable_number ("sigma_11");
302  sigma_vars[4] = stress_system.variable_number ("sigma_12");
303  sigma_vars[5] = stress_system.variable_number ("sigma_22");
304  unsigned int vonMises_var = stress_system.variable_number ("vonMises");
305 
306  // Storage for the stress dof indices on each element
307  std::vector<std::vector<dof_id_type>> dof_indices_var(system.n_vars());
308  std::vector<dof_id_type> stress_dof_indices_var;
309  std::vector<dof_id_type> vonmises_dof_indices_var;
310 
311  for (const auto & elem : mesh.active_local_element_ptr_range())
312  {
313  for (unsigned int var=0; var<3; var++)
314  dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
315 
316  const unsigned int n_var_dofs = dof_indices_var[0].size();
317 
318  fe->reinit (elem);
319 
320  std::vector<DenseMatrix<Number>> stress_tensor_qp(qrule.n_points());
321  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
322  {
323  stress_tensor_qp[qp].resize(3,3);
324 
325  // Row is variable u1, u2, or u3, column is x, y, or z
326  DenseMatrix<Number> grad_u(3,3);
327  for (unsigned int var_i=0; var_i<3; var_i++)
328  for (unsigned int var_j=0; var_j<3; var_j++)
329  for (unsigned int j=0; j<n_var_dofs; j++)
330  grad_u(var_i,var_j) += dphi[j][qp](var_j) * system.current_solution(dof_indices_var[var_i][j]);
331 
332  for (unsigned int var_i=0; var_i<3; var_i++)
333  for (unsigned int var_j=0; var_j<3; var_j++)
334  for (unsigned int k=0; k<3; k++)
335  for (unsigned int l=0; l<3; l++)
336  stress_tensor_qp[qp](var_i,var_j) += elasticity_tensor(var_i,var_j,k,l) * grad_u(k,l);
337  }
338 
339  stress_dof_map.dof_indices (elem, vonmises_dof_indices_var, vonMises_var);
340  std::vector<DenseMatrix<Number>> elem_sigma_vec(vonmises_dof_indices_var.size());
341  for (std::size_t index=0; index<elem_sigma_vec.size(); index++)
342  elem_sigma_vec[index].resize(3,3);
343 
344  // Below we project each component of the stress tensor onto a L2_LAGRANGE discretization.
345  // Note that this gives a discontinuous stress plot on element boundaries, which is
346  // appropriate. We then also get the von Mises stress from the projected stress tensor.
347  unsigned int stress_var_index = 0;
348  for (unsigned int var_i=0; var_i<3; var_i++)
349  for (unsigned int var_j=var_i; var_j<3; var_j++)
350  {
351  stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[stress_var_index]);
352 
353  const unsigned int n_proj_dofs = stress_dof_indices_var.size();
354 
355  DenseMatrix<Real> Me(n_proj_dofs, n_proj_dofs);
356  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
357  {
358  for(unsigned int i=0; i<n_proj_dofs; i++)
359  for(unsigned int j=0; j<n_proj_dofs; j++)
360  {
361  Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
362  }
363  }
364 
365  DenseVector<Number> Fe(n_proj_dofs);
366  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
367  for(unsigned int i=0; i<n_proj_dofs; i++)
368  {
369  Fe(i) += JxW[qp] * stress_tensor_qp[qp](var_i,var_j) * phi[i][qp];
370  }
371 
372  DenseVector<Number> projected_data;
373  Me.cholesky_solve(Fe, projected_data);
374 
375  for(unsigned int index=0; index<n_proj_dofs; index++)
376  {
377  dof_id_type dof_index = stress_dof_indices_var[index];
378  if ((stress_system.solution->first_local_index() <= dof_index) &&
379  (dof_index < stress_system.solution->last_local_index()))
380  stress_system.solution->set(dof_index, projected_data(index));
381 
382  elem_sigma_vec[index](var_i,var_j) = projected_data(index);
383  }
384 
385  stress_var_index++;
386  }
387 
388  for (std::size_t index=0; index<elem_sigma_vec.size(); index++)
389  {
390  elem_sigma_vec[index](1,0) = elem_sigma_vec[index](0,1);
391  elem_sigma_vec[index](2,0) = elem_sigma_vec[index](0,2);
392  elem_sigma_vec[index](2,1) = elem_sigma_vec[index](1,2);
393 
394  // Get the von Mises stress from the projected stress tensor
395  Number vonMises_value = std::sqrt(0.5*(pow(elem_sigma_vec[index](0,0) - elem_sigma_vec[index](1,1), 2.) +
396  pow(elem_sigma_vec[index](1,1) - elem_sigma_vec[index](2,2), 2.) +
397  pow(elem_sigma_vec[index](2,2) - elem_sigma_vec[index](0,0), 2.) +
398  6.*(pow(elem_sigma_vec[index](0,1), 2.) +
399  pow(elem_sigma_vec[index](1,2), 2.) +
400  pow(elem_sigma_vec[index](2,0), 2.))));
401 
402  dof_id_type dof_index = vonmises_dof_indices_var[index];
403 
404  if ((stress_system.solution->first_local_index() <= dof_index) &&
405  (dof_index < stress_system.solution->last_local_index()))
406  stress_system.solution->set(dof_index, vonMises_value);
407  }
408  }
409 
410  // Should call close and update when we set vector entries directly
411  stress_system.solution->close();
412  stress_system.update();
413  }
414 };
415 
416 
417 // Begin the main program.
418 int main (int argc, char ** argv)
419 {
420  // Initialize libMesh and any dependent libraries
421  LibMeshInit init (argc, argv);
422 
423  // This example requires a linear solver package.
424  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
425  "--enable-petsc, --enable-trilinos, or --enable-eigen");
426 
427  // Initialize the cantilever mesh
428  const unsigned int dim = 3;
429 
430  // Make sure libMesh was compiled for 3D
431  libmesh_example_requires(dim == LIBMESH_DIM, "3D support");
432 
433  // Create a 3D mesh distributed across the default MPI communicator.
434  Mesh mesh(init.comm(), dim);
436  32,
437  8,
438  4,
439  0., 1.*x_scaling,
440  0., 0.3,
441  0., 0.1,
442  HEX8);
443 
444  // Print information about the mesh to the screen.
445  mesh.print_info();
446 
447  // Let's add some node and edge boundary conditions
448  // Each processor should know about each boundary condition it can
449  // see, so we loop over all elements, not just local elements.
450  for (const auto & elem : mesh.element_ptr_range())
451  {
452  unsigned int
453  side_max_x = 0, side_min_y = 0,
454  side_max_y = 0, side_max_z = 0;
455 
456  bool
457  found_side_max_x = false, found_side_max_y = false,
458  found_side_min_y = false, found_side_max_z = false;
459 
460  for (auto side : elem->side_index_range())
461  {
462  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
463  {
464  side_max_x = side;
465  found_side_max_x = true;
466  }
467 
468  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MIN_Y))
469  {
470  side_min_y = side;
471  found_side_min_y = true;
472  }
473 
474  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Y))
475  {
476  side_max_y = side;
477  found_side_max_y = true;
478  }
479 
480  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_Z))
481  {
482  side_max_z = side;
483  found_side_max_z = true;
484  }
485  }
486 
487  // If elem has sides on boundaries
488  // BOUNDARY_ID_MAX_X, BOUNDARY_ID_MAX_Y, BOUNDARY_ID_MAX_Z
489  // then let's set a node boundary condition
490  if (found_side_max_x && found_side_max_y && found_side_max_z)
491  for (auto n : elem->node_index_range())
492  if (elem->is_node_on_side(n, side_max_x) &&
493  elem->is_node_on_side(n, side_max_y) &&
494  elem->is_node_on_side(n, side_max_z))
495  mesh.get_boundary_info().add_node(elem->node_ptr(n), NODE_BOUNDARY_ID);
496 
497 
498  // If elem has sides on boundaries
499  // BOUNDARY_ID_MAX_X and BOUNDARY_ID_MIN_Y
500  // then let's set an edge boundary condition
501  if (found_side_max_x && found_side_min_y)
502  for (auto e : elem->edge_index_range())
503  if (elem->is_edge_on_side(e, side_max_x) &&
504  elem->is_edge_on_side(e, side_min_y))
505  mesh.get_boundary_info().add_edge(elem, e, EDGE_BOUNDARY_ID);
506  }
507 
508  // Create an equation systems object.
509  EquationSystems equation_systems (mesh);
510 
511  // Declare the system and its variables.
512  // Create a system named "Elasticity"
513  LinearImplicitSystem & system =
514  equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
515 
516 #ifdef LIBMESH_HAVE_PETSC
517  // Attach a SolverConfiguration object to system.linear_solver
518  PetscLinearSolver<Number> * petsc_linear_solver =
519  cast_ptr<PetscLinearSolver<Number>*>(system.get_linear_solver());
520  libmesh_assert(petsc_linear_solver);
521  PetscSolverConfiguration petsc_solver_config(*petsc_linear_solver);
522  petsc_linear_solver->set_solver_configuration(petsc_solver_config);
523 #endif
524 
525  // Add three displacement variables, u and v, to the system
526  unsigned int u_var = system.add_variable("u", FIRST, LAGRANGE);
527  unsigned int v_var = system.add_variable("v", FIRST, LAGRANGE);
528  unsigned int w_var = system.add_variable("w", FIRST, LAGRANGE);
529 
530  LinearElasticity le(equation_systems);
531  system.attach_assemble_object(le);
532 
533  std::set<boundary_id_type> boundary_ids;
534  boundary_ids.insert(BOUNDARY_ID_MIN_X);
535  boundary_ids.insert(NODE_BOUNDARY_ID);
536  boundary_ids.insert(EDGE_BOUNDARY_ID);
537 
538  // Create a vector storing the variable numbers which the BC applies to
539  std::vector<unsigned int> variables;
540  variables.push_back(u_var);
541  variables.push_back(v_var);
542  variables.push_back(w_var);
543 
544  // Create a ZeroFunction to initialize dirichlet_bc
545  ZeroFunction<> zf;
546 
547  // Most DirichletBoundary users will want to supply a "locally
548  // indexed" functor
549  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
551 
552  // We must add the Dirichlet boundary condition _before_
553  // we call equation_systems.init()
554  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
555 
556  // Also, initialize an ExplicitSystem to store stresses
557  ExplicitSystem & stress_system =
558  equation_systems.add_system<ExplicitSystem> ("StressSystem");
559 
560  stress_system.add_variable("sigma_00", FIRST, L2_LAGRANGE);
561  stress_system.add_variable("sigma_01", FIRST, L2_LAGRANGE);
562  stress_system.add_variable("sigma_02", FIRST, L2_LAGRANGE);
563  stress_system.add_variable("sigma_11", FIRST, L2_LAGRANGE);
564  stress_system.add_variable("sigma_12", FIRST, L2_LAGRANGE);
565  stress_system.add_variable("sigma_22", FIRST, L2_LAGRANGE);
566  stress_system.add_variable("vonMises", FIRST, L2_LAGRANGE);
567 
568  // Initialize the data structures for the equation system.
569  equation_systems.init();
570 
571  // Print information about the system to the screen.
572  equation_systems.print_info();
573 
574  // Solve the system
575  system.solve();
576 
577  // Post-process the solution to compute the stresses
578  le.compute_stresses();
579 
580  // Plot the solution
581 #ifdef LIBMESH_HAVE_EXODUS_API
582 
583  // Use single precision in this case (reduces the size of the exodus file)
584  ExodusII_IO exo_io(mesh, /*single_precision=*/true);
585  exo_io.write_discontinuous_exodusII("displacement_and_stress.exo", equation_systems);
586 
587 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
588 
589  // All done.
590  return 0;
591 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
virtual void configure_solver()
Apply solver options to a particular solver.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
This is the EquationSystems class.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
void assemble()
Assemble the system matrix and right-hand side vector.
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
This class provides a specific system class.
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
unsigned short int side
Definition: xdr_io.C:49
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Definition: assembly.C:38
LinearElasticity(EquationSystems &es_in)
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
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.
Abstract base class to be used for system assembly.
Definition: system.h:119
int main(int argc, char **argv)
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
SolverPackage default_solver_package()
Definition: libmesh.C:995
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
unsigned short int variable_number(const std::string &var) const
Definition: system.C:1263
void add_node(const Node *node, const boundary_id_type id)
Add Node node with boundary id id to the boundary information data structures.
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.
virtual LinearSolver< Number > * get_linear_solver() const libmesh_override
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual SimpleRange< element_iterator > element_ptr_range()=0
PetscSolverConfiguration(PetscLinearSolver< Number > &petsc_linear_solver)
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
PetscLinearSolver< Number > & _petsc_linear_solver
Order default_quadrature_order() const
Definition: fe_type.h:332
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.
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
This class stores solver configuration data, e.g.
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
double pow(double a, int b)
Real elasticity_tensor(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 attach_assemble_object(Assembly &assemble)
Register a user object to use in assembling the system matrix and RHS.
Definition: system.C:1814
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
PetscErrorCode ierr
SparseMatrix< Number > * matrix
The system matrix.
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
const 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.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
Real kronecker_delta(unsigned int i, unsigned int j)
Kronecker delta function.
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
Definition: dense_matrix.C:911
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 write_discontinuous_exodusII(const std::string &name, const EquationSystems &es, const std::set< std::string > *system_names=libmesh_nullptr)
Writes a exodusII file with discontinuous data.
Definition: exodusII_io.C:89
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 add_edge(const dof_id_type elem, const unsigned short int edge, const boundary_id_type id)
Add edge edge of element number elem with boundary id id to the boundary information data structure...
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.