1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 #include "libmesh/enum_solver_package.h"
60
61 // Bring in everything from the libMesh namespace
62 using namespace libMesh;
63
64 // Matrix and right-hand side assemble
65 void assemble_elasticity(EquationSystems & es,
66 const std::string & system_name);
67
68 // Define the elasticity tensor, which is a fourth-order tensor
69 // i.e. it has four indices i, j, k, l
70 Real eval_elasticity_tensor(unsigned int i,
71 unsigned int j,
72 unsigned int k,
73 unsigned int l);
74
75 // Begin the main program.
main(int argc,char ** argv)76 int main (int argc, char ** argv)
77 {
78 // Initialize libMesh and any dependent libraries
79 LibMeshInit init (argc, argv);
80
81 // This example requires a linear solver package.
82 libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
83 "--enable-petsc, --enable-trilinos, or --enable-eigen");
84
85 // Initialize the cantilever mesh
86 const unsigned int dim = 2;
87
88 // Skip this 2D example if libMesh was compiled as 1D-only.
89 libmesh_example_requires(dim <= LIBMESH_DIM, "2D support");
90
91 // We use Dirichlet boundary conditions here
92 #ifndef LIBMESH_ENABLE_DIRICHLET
93 libmesh_example_requires(false, "--enable-dirichlet");
94 #endif
95
96 // Create a 2D mesh distributed across the default MPI communicator.
97 Mesh mesh(init.comm(), dim);
98 MeshTools::Generation::build_square (mesh,
99 50, 10,
100 0., 1.,
101 0., 0.2,
102 QUAD9);
103
104
105 // Print information about the mesh to the screen.
106 mesh.print_info();
107
108 // Create an equation systems object.
109 EquationSystems equation_systems (mesh);
110
111 // Declare the system and its variables.
112 // Create a system named "Elasticity"
113 LinearImplicitSystem & system =
114 equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
115
116 system.attach_assemble_function (assemble_elasticity);
117
118 #ifdef LIBMESH_ENABLE_DIRICHLET
119 // Add two displacement variables, u and v, to the system
120 unsigned int u_var = system.add_variable("u", SECOND, LAGRANGE);
121 unsigned int v_var = system.add_variable("v", SECOND, LAGRANGE);
122
123 // Construct a Dirichlet boundary condition object
124 // We impose a "clamped" boundary condition on the
125 // "left" boundary, i.e. bc_id = 3
126 std::set<boundary_id_type> boundary_ids;
127 boundary_ids.insert(3);
128
129 // Create a vector storing the variable numbers which the BC applies to
130 std::vector<unsigned int> variables(2);
131 variables[0] = u_var; variables[1] = v_var;
132
133 // Create a ZeroFunction to initialize dirichlet_bc
134 ZeroFunction<> zf;
135
136 // Most DirichletBoundary users will want to supply a "locally
137 // indexed" functor
138 DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
139 LOCAL_VARIABLE_ORDER);
140
141 // We must add the Dirichlet boundary condition _before_
142 // we call equation_systems.init()
143 system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
144 #endif // LIBMESH_ENABLE_DIRICHLET
145
146 // Initialize the data structures for the equation system.
147 equation_systems.init();
148
149 // Print information about the system to the screen.
150 equation_systems.print_info();
151
152 // Solve the system
153 system.solve();
154
155 // Plot the solution
156 #ifdef LIBMESH_HAVE_EXODUS_API
157 ExodusII_IO (mesh).write_equation_systems("displacement.e", equation_systems);
158 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
159
160 // All done.
161 return 0;
162 }
163
164
assemble_elasticity(EquationSystems & es,const std::string & libmesh_dbg_var (system_name))165 void assemble_elasticity(EquationSystems & es,
166 const std::string & libmesh_dbg_var(system_name))
167 {
168 libmesh_assert_equal_to (system_name, "Elasticity");
169
170 const MeshBase & mesh = es.get_mesh();
171
172 const unsigned int dim = mesh.mesh_dimension();
173
174 LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Elasticity");
175
176 const unsigned int u_var = system.variable_number ("u");
177 const unsigned int v_var = system.variable_number ("v");
178
179 const DofMap & dof_map = system.get_dof_map();
180 FEType fe_type = dof_map.variable_type(0);
181 std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
182 QGauss qrule (dim, fe_type.default_quadrature_order());
183 fe->attach_quadrature_rule (&qrule);
184
185 std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
186 QGauss qface(dim-1, fe_type.default_quadrature_order());
187 fe_face->attach_quadrature_rule (&qface);
188
189 const std::vector<Real> & JxW = fe->get_JxW();
190 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
191
192 DenseMatrix<Number> Ke;
193 DenseVector<Number> Fe;
194
195 DenseSubMatrix<Number>
196 Kuu(Ke), Kuv(Ke),
197 Kvu(Ke), Kvv(Ke);
198
199 DenseSubVector<Number>
200 Fu(Fe),
201 Fv(Fe);
202
203 std::vector<dof_id_type> dof_indices;
204 std::vector<dof_id_type> dof_indices_u;
205 std::vector<dof_id_type> dof_indices_v;
206
207 for (const auto & elem : mesh.active_local_element_ptr_range())
208 {
209 dof_map.dof_indices (elem, dof_indices);
210 dof_map.dof_indices (elem, dof_indices_u, u_var);
211 dof_map.dof_indices (elem, dof_indices_v, v_var);
212
213 const unsigned int n_dofs = dof_indices.size();
214 const unsigned int n_u_dofs = dof_indices_u.size();
215 const unsigned int n_v_dofs = dof_indices_v.size();
216
217 fe->reinit (elem);
218
219 Ke.resize (n_dofs, n_dofs);
220 Fe.resize (n_dofs);
221
222 Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
223 Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
224
225 Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
226 Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
227
228 Fu.reposition (u_var*n_u_dofs, n_u_dofs);
229 Fv.reposition (v_var*n_u_dofs, n_v_dofs);
230
231 for (unsigned int qp=0; qp<qrule.n_points(); qp++)
232 {
233 for (unsigned int i=0; i<n_u_dofs; i++)
234 for (unsigned int j=0; j<n_u_dofs; j++)
235 {
236 // Tensor indices
237 unsigned int C_i, C_j, C_k, C_l;
238 C_i=0, C_k=0;
239
240 C_j=0, C_l=0;
241 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));
242
243 C_j=1, C_l=0;
244 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));
245
246 C_j=0, C_l=1;
247 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));
248
249 C_j=1, C_l=1;
250 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));
251 }
252
253 for (unsigned int i=0; i<n_u_dofs; i++)
254 for (unsigned int j=0; j<n_v_dofs; j++)
255 {
256 // Tensor indices
257 unsigned int C_i, C_j, C_k, C_l;
258 C_i=0, C_k=1;
259
260 C_j=0, C_l=0;
261 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));
262
263 C_j=1, C_l=0;
264 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));
265
266 C_j=0, C_l=1;
267 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));
268
269 C_j=1, C_l=1;
270 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));
271 }
272
273 for (unsigned int i=0; i<n_v_dofs; i++)
274 for (unsigned int j=0; j<n_u_dofs; j++)
275 {
276 // Tensor indices
277 unsigned int C_i, C_j, C_k, C_l;
278 C_i=1, C_k=0;
279
280 C_j=0, C_l=0;
281 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));
282
283 C_j=1, C_l=0;
284 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));
285
286 C_j=0, C_l=1;
287 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));
288
289 C_j=1, C_l=1;
290 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));
291 }
292
293 for (unsigned int i=0; i<n_v_dofs; i++)
294 for (unsigned int j=0; j<n_v_dofs; j++)
295 {
296 // Tensor indices
297 unsigned int C_i, C_j, C_k, C_l;
298 C_i=1, C_k=1;
299
300 C_j=0, C_l=0;
301 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));
302
303 C_j=1, C_l=0;
304 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));
305
306 C_j=0, C_l=1;
307 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));
308
309 C_j=1, C_l=1;
310 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));
311 }
312 }
313
314 {
315 for (auto side : elem->side_index_range())
316 if (elem->neighbor_ptr(side) == nullptr)
317 {
318 const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
319 const std::vector<Real> & JxW_face = fe_face->get_JxW();
320
321 fe_face->reinit(elem, side);
322
323 if (mesh.get_boundary_info().has_boundary_id (elem, side, 1)) // Apply a traction on the right side
324 {
325 for (unsigned int qp=0; qp<qface.n_points(); qp++)
326 for (unsigned int i=0; i<n_v_dofs; i++)
327 Fv(i) += JxW_face[qp] * (-1.) * phi_face[i][qp];
328 }
329 }
330 }
331
332 dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
333
334 system.matrix->add_matrix (Ke, dof_indices);
335 system.rhs->add_vector (Fe, dof_indices);
336 }
337 }
338
eval_elasticity_tensor(unsigned int i,unsigned int j,unsigned int k,unsigned int l)339 Real eval_elasticity_tensor(unsigned int i,
340 unsigned int j,
341 unsigned int k,
342 unsigned int l)
343 {
344 // Define the Poisson ratio
345 const Real nu = 0.3;
346
347 // Define the Lame constants (lambda_1 and lambda_2) based on Poisson ratio
348 const Real lambda_1 = nu / ((1. + nu) * (1. - 2.*nu));
349 const Real lambda_2 = 0.5 / (1 + nu);
350
351 // Define the Kronecker delta functions that we need here
352 Real delta_ij = (i == j) ? 1. : 0.;
353 Real delta_il = (i == l) ? 1. : 0.;
354 Real delta_ik = (i == k) ? 1. : 0.;
355 Real delta_jl = (j == l) ? 1. : 0.;
356 Real delta_jk = (j == k) ? 1. : 0.;
357 Real delta_kl = (k == l) ? 1. : 0.;
358
359 return lambda_1 * delta_ij * delta_kl + lambda_2 * (delta_ik * delta_jl + delta_il * delta_jk);
360 }
361