libMesh
systems_of_equations_ex5.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
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 5 - Linear Elastic Cantilever with Constraint </h1>
21 // \author David Knezevic
22 // \date 2012
23 //
24 // In this example we extend systems_of_equations_ex4 to enforce a constraint.
25 // We apply a uniform load on the top surface of the cantilever, and we
26 // determine the traction on the right boundary in order to obtain zero
27 // average vertical displacement on the right boundary of the domain.
28 //
29 // This constraint is enforced via a Lagrange multiplier (SCALAR variable).
30 // The system we solve, therefore, is of the form:
31 // a(u,v) + \lambda g(v) = f(v)
32 // g(u) = 0
33 // Here \lambda tells us the traction required to satisfy the constraint.
34
35 // C++ include files that we need
36 #include <iostream>
37 #include <algorithm>
38 #include <math.h>
39
40 // libMesh includes
41 #include "libmesh/libmesh.h"
42 #include "libmesh/mesh.h"
43 #include "libmesh/mesh_generation.h"
44 #include "libmesh/exodusII_io.h"
45 #include "libmesh/linear_implicit_system.h"
46 #include "libmesh/equation_systems.h"
47 #include "libmesh/fe.h"
49 #include "libmesh/dof_map.h"
50 #include "libmesh/sparse_matrix.h"
51 #include "libmesh/numeric_vector.h"
52 #include "libmesh/dense_matrix.h"
53 #include "libmesh/dense_submatrix.h"
54 #include "libmesh/dense_vector.h"
55 #include "libmesh/dense_subvector.h"
56 #include "libmesh/perf_log.h"
57 #include "libmesh/elem.h"
58 #include "libmesh/boundary_info.h"
59 #include "libmesh/zero_function.h"
60 #include "libmesh/dirichlet_boundaries.h"
61 #include "libmesh/string_to_enum.h"
62 #include "libmesh/getpot.h"
63
64 // Bring in everything from the libMesh namespace
65 using namespace libMesh;
66
67 // Matrix and right-hand side assemble
69  const std::string & system_name);
70
71 // Define the elasticity tensor, which is a fourth-order tensor
72 // i.e. it has four indices i, j, k, l
73 Real eval_elasticity_tensor(unsigned int i,
74  unsigned int j,
75  unsigned int k,
76  unsigned int l);
77
78 // Begin the main program.
79 int main (int argc, char ** argv)
80 {
81  // Initialize libMesh and any dependent libraries
82  LibMeshInit init (argc, argv);
83
84  // This example requires a linear solver package.
85  libmesh_example_requires(libMesh::default_solver_package() != INVALID_SOLVER_PACKAGE,
86  "--enable-petsc, --enable-trilinos, or --enable-eigen");
87
88  // This example NaNs with the Eigen sparse linear solvers
89  libmesh_example_requires(libMesh::default_solver_package() != EIGEN_SOLVERS, "--enable-petsc or --enable-laspack");
90
91  // Initialize the cantilever mesh
92  const unsigned int dim = 2;
93
94  // Skip this 2D example if libMesh was compiled as 1D-only.
95  libmesh_example_requires(dim <= LIBMESH_DIM, "2D support");
96
97  // Create a 2D mesh distributed across the default MPI communicator.
98  Mesh mesh(init.comm(), dim);
100  50, 10,
101  0., 1.,
102  0., 0.2,
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 =
115
116  // Add two displacement variables, u and v, to the system
117  unsigned int u_var = system.add_variable("u", SECOND, LAGRANGE);
118  unsigned int v_var = system.add_variable("v", SECOND, LAGRANGE);
119
120  // Add a SCALAR variable for the Lagrange multiplier to enforce our constraint
122
124
125  // Construct a Dirichlet boundary condition object
126  // We impose a "clamped" boundary condition on the
127  // "left" boundary, i.e. bc_id = 3
128  std::set<boundary_id_type> boundary_ids;
129  boundary_ids.insert(3);
130
131  // Create a vector storing the variable numbers which the BC applies to
132  std::vector<unsigned int> variables(2);
133  variables[0] = u_var;
134  variables[1] = v_var;
135
136  // Create a ZeroFunction to initialize dirichlet_bc
137  ZeroFunction<> zf;
138
139  // Most DirichletBoundary users will want to supply a "locally
140  // indexed" functor
141  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
143
144  // We must add the Dirichlet boundary condition _before_
145  // we call equation_systems.init()
147
148  // Initialize the data structures for the equation system.
149  equation_systems.init();
150
151  // Print information about the system to the screen.
152  equation_systems.print_info();
153
154  // Solve the system
155  system.solve();
156
157  // Plot the solution
158 #ifdef LIBMESH_HAVE_EXODUS_API
159  ExodusII_IO (mesh).write_equation_systems("displacement.e", equation_systems);
160 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
161
162  // All done.
163  return 0;
164 }
165
166
168  const std::string & libmesh_dbg_var(system_name))
169 {
170  libmesh_assert_equal_to (system_name, "Elasticity");
171
172  const MeshBase & mesh = es.get_mesh();
173
174  const unsigned int dim = mesh.mesh_dimension();
175
176  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Elasticity");
177
178  const unsigned int u_var = system.variable_number ("u");
179  const unsigned int v_var = system.variable_number ("v");
180  const unsigned int lambda_var = system.variable_number ("lambda");
181
182  const DofMap & dof_map = system.get_dof_map();
183  FEType fe_type = dof_map.variable_type(0);
184  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
187
188  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
191
192  const std::vector<Real> & JxW = fe->get_JxW();
193  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
194
197
199  Kuu(Ke), Kuv(Ke),
200  Kvu(Ke), Kvv(Ke);
201  DenseSubMatrix<Number> Klambda_v(Ke), Kv_lambda(Ke);
202
204  Fu(Fe),
205  Fv(Fe);
206
207  std::vector<dof_id_type> dof_indices;
208  std::vector<dof_id_type> dof_indices_u;
209  std::vector<dof_id_type> dof_indices_v;
210  std::vector<dof_id_type> dof_indices_lambda;
211
212  for (const auto & elem : mesh.active_local_element_ptr_range())
213  {
214  dof_map.dof_indices (elem, dof_indices);
215  dof_map.dof_indices (elem, dof_indices_u, u_var);
216  dof_map.dof_indices (elem, dof_indices_v, v_var);
217  dof_map.dof_indices (elem, dof_indices_lambda, lambda_var);
218
219  const unsigned int n_dofs = dof_indices.size();
220  const unsigned int n_u_dofs = dof_indices_u.size();
221  const unsigned int n_v_dofs = dof_indices_v.size();
222  const unsigned int n_lambda_dofs = dof_indices_lambda.size();
223
224  fe->reinit (elem);
225
226  Ke.resize (n_dofs, n_dofs);
227  Fe.resize (n_dofs);
228
229  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
230  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
231
232  Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
233  Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
234
235  // Also, add a row and a column to enforce the constraint
236  Kv_lambda.reposition (v_var*n_u_dofs, v_var*n_u_dofs+n_v_dofs, n_v_dofs, 1);
237  Klambda_v.reposition (v_var*n_v_dofs+n_v_dofs, v_var*n_v_dofs, 1, n_v_dofs);
238
239  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
240  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
241
242  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
243  {
244  for (unsigned int i=0; i<n_u_dofs; i++)
245  for (unsigned int j=0; j<n_u_dofs; j++)
246  {
247  // Tensor indices
248  unsigned int C_i, C_j, C_k, C_l;
249  C_i=0, C_k=0;
250
251  C_j=0, C_l=0;
252  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));
253
254  C_j=1, C_l=0;
255  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));
256
257  C_j=0, C_l=1;
258  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));
259
260  C_j=1, C_l=1;
261  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));
262  }
263
264  for (unsigned int i=0; i<n_u_dofs; i++)
265  for (unsigned int j=0; j<n_v_dofs; j++)
266  {
267  // Tensor indices
268  unsigned int C_i, C_j, C_k, C_l;
269  C_i=0, C_k=1;
270
271  C_j=0, C_l=0;
272  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));
273
274  C_j=1, C_l=0;
275  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));
276
277  C_j=0, C_l=1;
278  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));
279
280  C_j=1, C_l=1;
281  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));
282  }
283
284  for (unsigned int i=0; i<n_v_dofs; i++)
285  for (unsigned int j=0; j<n_u_dofs; j++)
286  {
287  // Tensor indices
288  unsigned int C_i, C_j, C_k, C_l;
289  C_i=1, C_k=0;
290
291  C_j=0, C_l=0;
292  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));
293
294  C_j=1, C_l=0;
295  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));
296
297  C_j=0, C_l=1;
298  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));
299
300  C_j=1, C_l=1;
301  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));
302  }
303
304  for (unsigned int i=0; i<n_v_dofs; i++)
305  for (unsigned int j=0; j<n_v_dofs; j++)
306  {
307  // Tensor indices
308  unsigned int C_i, C_j, C_k, C_l;
309  C_i=1, C_k=1;
310
311  C_j=0, C_l=0;
312  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));
313
314  C_j=1, C_l=0;
315  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));
316
317  C_j=0, C_l=1;
318  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));
319
320  C_j=1, C_l=1;
321  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));
322  }
323  }
324
325  {
326  std::vector<boundary_id_type> bc_ids;
327  for (auto side : elem->side_index_range())
328  if (elem->neighbor_ptr(side) == libmesh_nullptr)
329  {
330  mesh.get_boundary_info().boundary_ids (elem, side, bc_ids);
331
332  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
333  const std::vector<Real> & JxW_face = fe_face->get_JxW();
334
335  fe_face->reinit(elem, side);
336
337  for (std::vector<boundary_id_type>::const_iterator b =
338  bc_ids.begin(); b != bc_ids.end(); ++b)
339  {
340  const boundary_id_type bc_id = *b;
341  for (unsigned int qp=0; qp<qface.n_points(); qp++)
342  {
344  if (bc_id == 2)
345  for (unsigned int i=0; i<n_v_dofs; i++)
346  Fv(i) += JxW_face[qp] * (-1.) * phi_face[i][qp];
347
348  // Add the constraint contributions
349  if (bc_id == 1)
350  {
351  for (unsigned int i=0; i<n_v_dofs; i++)
352  for (unsigned int j=0; j<n_lambda_dofs; j++)
353  Kv_lambda(i,j) += JxW_face[qp] * (-1.) * phi_face[i][qp];
354
355  for (unsigned int i=0; i<n_lambda_dofs; i++)
356  for (unsigned int j=0; j<n_v_dofs; j++)
357  Klambda_v(i,j) += JxW_face[qp] * (-1.) * phi_face[j][qp];
358  }
359  }
360  }
361  }
362  }
363
364  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
365
368  }
369 }
370
372  unsigned int j,
373  unsigned int k,
374  unsigned int l)
375 {
376  // Define the Poisson ratio
377  const Real nu = 0.3;
378
379  // Define the Lame constants (lambda_1 and lambda_2) based on Poisson ratio
380  const Real lambda_1 = nu / ((1. + nu) * (1. - 2.*nu));
381  const Real lambda_2 = 0.5 / (1 + nu);
382
383  // Define the Kronecker delta functions that we need here
384  Real delta_ij = (i == j) ? 1. : 0.;
385  Real delta_il = (i == l) ? 1. : 0.;
386  Real delta_ik = (i == k) ? 1. : 0.;
387  Real delta_jl = (j == l) ? 1. : 0.;
388  Real delta_jk = (j == k) ? 1. : 0.;
389  Real delta_kl = (k == l) ? 1. : 0.;
390
391  return lambda_1 * delta_ij * delta_kl + lambda_2 * (delta_ik * delta_jl + delta_il * delta_jk);
392 }
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.
boundary_id_type bc_id
Definition: xdr_io.C:50
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
std::vector< boundary_id_type > boundary_ids(const Node *node) const
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.
const DofMap & get_dof_map() const
Definition: system.h:2030
int8_t boundary_id_type
Definition: id_types.h:51
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
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.
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
Adds a copy of the specified Dirichlet boundary to the system.
Real eval_elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
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
void assemble_elasticity(EquationSystems &es, const std::string &system_name)
int main(int argc, char **argv)
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.