libMesh
systems_of_equations_ex4.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 4 - Linear Elastic Cantilever </h1>
21 // \author David Knezevic
22 // \date 2012
23 //
24 // In this example we model a homogeneous isotropic cantilever
25 // using the equations of linear elasticity. We set the Poisson ratio to
26 // \nu = 0.3 and clamp the left boundary and apply a vertical load at the
27 // right boundary.
28 
29 
30 // C++ include files that we need
31 #include <iostream>
32 #include <algorithm>
33 #include <math.h>
34 
35 // libMesh includes
36 #include "libmesh/libmesh.h"
37 #include "libmesh/mesh.h"
38 #include "libmesh/mesh_generation.h"
39 #include "libmesh/exodusII_io.h"
40 #include "libmesh/gnuplot_io.h"
41 #include "libmesh/linear_implicit_system.h"
42 #include "libmesh/equation_systems.h"
43 #include "libmesh/fe.h"
44 #include "libmesh/quadrature_gauss.h"
45 #include "libmesh/dof_map.h"
46 #include "libmesh/sparse_matrix.h"
47 #include "libmesh/numeric_vector.h"
48 #include "libmesh/dense_matrix.h"
49 #include "libmesh/dense_submatrix.h"
50 #include "libmesh/dense_vector.h"
51 #include "libmesh/dense_subvector.h"
52 #include "libmesh/perf_log.h"
53 #include "libmesh/elem.h"
54 #include "libmesh/boundary_info.h"
55 #include "libmesh/zero_function.h"
56 #include "libmesh/dirichlet_boundaries.h"
57 #include "libmesh/string_to_enum.h"
58 #include "libmesh/getpot.h"
59 
60 // Bring in everything from the libMesh namespace
61 using namespace libMesh;
62 
63 // Matrix and right-hand side assemble
65  const std::string & system_name);
66 
67 // Define the elasticity tensor, which is a fourth-order tensor
68 // i.e. it has four indices i, j, k, l
69 Real eval_elasticity_tensor(unsigned int i,
70  unsigned int j,
71  unsigned int k,
72  unsigned int l);
73 
74 // Begin the main program.
75 int main (int argc, char ** argv)
76 {
77  // Initialize libMesh and any dependent libraries
78  LibMeshInit init (argc, argv);
79 
80  // This example requires a linear solver package.
81  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
82  "--enable-petsc, --enable-trilinos, or --enable-eigen");
83 
84  // Initialize the cantilever mesh
85  const unsigned int dim = 2;
86 
87  // Skip this 2D example if libMesh was compiled as 1D-only.
88  libmesh_example_requires(dim <= LIBMESH_DIM, "2D support");
89 
90  // Create a 2D mesh distributed across the default MPI communicator.
91  Mesh mesh(init.comm(), dim);
93  50, 10,
94  0., 1.,
95  0., 0.2,
96  QUAD9);
97 
98 
99  // Print information about the mesh to the screen.
100  mesh.print_info();
101 
102  // Create an equation systems object.
103  EquationSystems equation_systems (mesh);
104 
105  // Declare the system and its variables.
106  // Create a system named "Elasticity"
107  LinearImplicitSystem & system =
108  equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
109 
110  // Add two displacement variables, u and v, to the system
111  unsigned int u_var = system.add_variable("u", SECOND, LAGRANGE);
112  unsigned int v_var = system.add_variable("v", SECOND, LAGRANGE);
113 
115 
116  // Construct a Dirichlet boundary condition object
117  // We impose a "clamped" boundary condition on the
118  // "left" boundary, i.e. bc_id = 3
119  std::set<boundary_id_type> boundary_ids;
120  boundary_ids.insert(3);
121 
122  // Create a vector storing the variable numbers which the BC applies to
123  std::vector<unsigned int> variables(2);
124  variables[0] = u_var; variables[1] = v_var;
125 
126  // Create a ZeroFunction to initialize dirichlet_bc
127  ZeroFunction<> zf;
128 
129  // Most DirichletBoundary users will want to supply a "locally
130  // indexed" functor
131  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
133 
134  // We must add the Dirichlet boundary condition _before_
135  // we call equation_systems.init()
136  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
137 
138  // Initialize the data structures for the equation system.
139  equation_systems.init();
140 
141  // Print information about the system to the screen.
142  equation_systems.print_info();
143 
144  // Solve the system
145  system.solve();
146 
147  // Plot the solution
148 #ifdef LIBMESH_HAVE_EXODUS_API
149  ExodusII_IO (mesh).write_equation_systems("displacement.e", equation_systems);
150 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
151 
152  // All done.
153  return 0;
154 }
155 
156 
158  const std::string & libmesh_dbg_var(system_name))
159 {
160  libmesh_assert_equal_to (system_name, "Elasticity");
161 
162  const MeshBase & mesh = es.get_mesh();
163 
164  const unsigned int dim = mesh.mesh_dimension();
165 
166  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Elasticity");
167 
168  const unsigned int u_var = system.variable_number ("u");
169  const unsigned int v_var = system.variable_number ("v");
170 
171  const DofMap & dof_map = system.get_dof_map();
172  FEType fe_type = dof_map.variable_type(0);
173  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
174  QGauss qrule (dim, fe_type.default_quadrature_order());
175  fe->attach_quadrature_rule (&qrule);
176 
177  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
178  QGauss qface(dim-1, fe_type.default_quadrature_order());
179  fe_face->attach_quadrature_rule (&qface);
180 
181  const std::vector<Real> & JxW = fe->get_JxW();
182  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
183 
186 
188  Kuu(Ke), Kuv(Ke),
189  Kvu(Ke), Kvv(Ke);
190 
192  Fu(Fe),
193  Fv(Fe);
194 
195  std::vector<dof_id_type> dof_indices;
196  std::vector<dof_id_type> dof_indices_u;
197  std::vector<dof_id_type> dof_indices_v;
198 
199  for (const auto & elem : mesh.active_local_element_ptr_range())
200  {
201  dof_map.dof_indices (elem, dof_indices);
202  dof_map.dof_indices (elem, dof_indices_u, u_var);
203  dof_map.dof_indices (elem, dof_indices_v, v_var);
204 
205  const unsigned int n_dofs = dof_indices.size();
206  const unsigned int n_u_dofs = dof_indices_u.size();
207  const unsigned int n_v_dofs = dof_indices_v.size();
208 
209  fe->reinit (elem);
210 
211  Ke.resize (n_dofs, n_dofs);
212  Fe.resize (n_dofs);
213 
214  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
215  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
216 
217  Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
218  Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
219 
220  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
221  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
222 
223  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
224  {
225  for (unsigned int i=0; i<n_u_dofs; i++)
226  for (unsigned int j=0; j<n_u_dofs; j++)
227  {
228  // Tensor indices
229  unsigned int C_i, C_j, C_k, C_l;
230  C_i=0, C_k=0;
231 
232  C_j=0, C_l=0;
233  Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
234 
235  C_j=1, C_l=0;
236  Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
237 
238  C_j=0, C_l=1;
239  Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
240 
241  C_j=1, C_l=1;
242  Kuu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
243  }
244 
245  for (unsigned int i=0; i<n_u_dofs; i++)
246  for (unsigned int j=0; j<n_v_dofs; j++)
247  {
248  // Tensor indices
249  unsigned int C_i, C_j, C_k, C_l;
250  C_i=0, C_k=1;
251 
252  C_j=0, C_l=0;
253  Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
254 
255  C_j=1, C_l=0;
256  Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
257 
258  C_j=0, C_l=1;
259  Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
260 
261  C_j=1, C_l=1;
262  Kuv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
263  }
264 
265  for (unsigned int i=0; i<n_v_dofs; i++)
266  for (unsigned int j=0; j<n_u_dofs; j++)
267  {
268  // Tensor indices
269  unsigned int C_i, C_j, C_k, C_l;
270  C_i=1, C_k=0;
271 
272  C_j=0, C_l=0;
273  Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
274 
275  C_j=1, C_l=0;
276  Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
277 
278  C_j=0, C_l=1;
279  Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
280 
281  C_j=1, C_l=1;
282  Kvu(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
283  }
284 
285  for (unsigned int i=0; i<n_v_dofs; i++)
286  for (unsigned int j=0; j<n_v_dofs; j++)
287  {
288  // Tensor indices
289  unsigned int C_i, C_j, C_k, C_l;
290  C_i=1, C_k=1;
291 
292  C_j=0, C_l=0;
293  Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
294 
295  C_j=1, C_l=0;
296  Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
297 
298  C_j=0, C_l=1;
299  Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
300 
301  C_j=1, C_l=1;
302  Kvv(i,j) += JxW[qp]*(eval_elasticity_tensor(C_i, C_j, C_k, C_l) * dphi[i][qp](C_j)*dphi[j][qp](C_l));
303  }
304  }
305 
306  {
307  for (auto side : elem->side_index_range())
308  if (elem->neighbor_ptr(side) == libmesh_nullptr)
309  {
310  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
311  const std::vector<Real> & JxW_face = fe_face->get_JxW();
312 
313  fe_face->reinit(elem, side);
314 
315  if (mesh.get_boundary_info().has_boundary_id (elem, side, 1)) // Apply a traction on the right side
316  {
317  for (unsigned int qp=0; qp<qface.n_points(); qp++)
318  for (unsigned int i=0; i<n_v_dofs; i++)
319  Fv(i) += JxW_face[qp] * (-1.) * phi_face[i][qp];
320  }
321  }
322  }
323 
324  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
325 
326  system.matrix->add_matrix (Ke, dof_indices);
327  system.rhs->add_vector (Fe, dof_indices);
328  }
329 }
330 
332  unsigned int j,
333  unsigned int k,
334  unsigned int l)
335 {
336  // Define the Poisson ratio
337  const Real nu = 0.3;
338 
339  // Define the Lame constants (lambda_1 and lambda_2) based on Poisson ratio
340  const Real lambda_1 = nu / ((1. + nu) * (1. - 2.*nu));
341  const Real lambda_2 = 0.5 / (1 + nu);
342 
343  // Define the Kronecker delta functions that we need here
344  Real delta_ij = (i == j) ? 1. : 0.;
345  Real delta_il = (i == l) ? 1. : 0.;
346  Real delta_ik = (i == k) ? 1. : 0.;
347  Real delta_jl = (j == l) ? 1. : 0.;
348  Real delta_jk = (j == k) ? 1. : 0.;
349  Real delta_kl = (k == l) ? 1. : 0.;
350 
351  return lambda_1 * delta_ij * delta_kl + lambda_2 * (delta_ik * delta_jl + delta_il * delta_jk);
352 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
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
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.
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
unsigned short int side
Definition: xdr_io.C:49
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...
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
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.
This is the MeshBase class.
Definition: mesh_base.h:68
Real eval_elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
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 reposition(const unsigned int ioff, const unsigned int n)
Changes the location of the subvector in the parent vector.
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.
int main(int argc, char **argv)
void assemble_elasticity(EquationSystems &es, const std::string &system_name)
const DofMap & get_dof_map() const
Definition: system.h:2030
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:1795
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
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.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
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
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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.
void reposition(const unsigned int ioff, const unsigned int joff, const unsigned int new_m, const unsigned int new_n)
Changes the location of the submatrix in the parent matrix.
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual void solve() libmesh_override
Assembles & solves the linear system A*x=b.
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
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.