libMesh
Functions
systems_of_equations_ex5.C File Reference

Go to the source code of this file.

Functions

void assemble_elasticity (EquationSystems &es, const std::string &system_name)
 
Real eval_elasticity_tensor (unsigned int i, unsigned int j, unsigned int k, unsigned int l)
 
int main (int argc, char **argv)
 
void assemble_elasticity (EquationSystems &es, const std::string &libmesh_dbg_var(system_name))
 

Function Documentation

void assemble_elasticity ( EquationSystems es,
const std::string &  system_name 
)

Referenced by main().

void assemble_elasticity ( EquationSystems es,
const std::string &  libmesh_dbg_varsystem_name 
)

Definition at line 167 of file systems_of_equations_ex5.C.

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), bc_id, libMesh::BoundaryInfo::boundary_ids(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DofMap::constrain_element_matrix_and_vector(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), eval_elasticity_tensor(), libMesh::MeshBase::get_boundary_info(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libmesh_nullptr, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::DenseSubVector< T >::reposition(), libMesh::DenseSubMatrix< T >::reposition(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, side, libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

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));
185  QGauss qrule (dim, fe_type.default_quadrature_order());
186  fe->attach_quadrature_rule (&qrule);
187 
188  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
189  QGauss qface(dim-1, fe_type.default_quadrature_order());
190  fe_face->attach_quadrature_rule (&qface);
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  {
343  // Add the loading
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 
366  system.matrix->add_matrix (Ke, dof_indices);
367  system.rhs->add_vector (Fe, dof_indices);
368  }
369 }
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
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
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 short int side
Definition: xdr_io.C:49
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
NumericVector< Number > * rhs
The system matrix.
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
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
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
Order default_quadrature_order() const
Definition: fe_type.h:332
Defines a dense submatrix for use in Finite Element-type computations.
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
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
Real eval_elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
const T_sys & get_system(const std::string &name) const
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
Real eval_elasticity_tensor ( unsigned int  i,
unsigned int  j,
unsigned int  k,
unsigned int  l 
)

Definition at line 371 of file systems_of_equations_ex5.C.

References libMesh::Real.

Referenced by assemble_elasticity().

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 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
int main ( int  argc,
char **  argv 
)

Definition at line 79 of file systems_of_equations_ex5.C.

References libMesh::DofMap::add_dirichlet_boundary(), libMesh::EquationSystems::add_system(), libMesh::System::add_variable(), assemble_elasticity(), libMesh::System::attach_assemble_function(), libMesh::MeshTools::Generation::build_square(), libMesh::LibMeshInit::comm(), libMesh::default_solver_package(), dim, libMesh::EIGEN_SOLVERS, libMesh::FIRST, libMesh::System::get_dof_map(), libMesh::TriangleWrapper::init(), libMesh::EquationSystems::init(), libMesh::INVALID_SOLVER_PACKAGE, libMesh::LAGRANGE, libMesh::LOCAL_VARIABLE_ORDER, mesh, libMesh::EquationSystems::print_info(), libMesh::MeshBase::print_info(), libMesh::QUAD9, libMesh::SCALAR, libMesh::SECOND, libMesh::LinearImplicitSystem::solve(), and libMesh::MeshOutput< MT >::write_equation_systems().

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,
103  QUAD9);
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  // 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
121  system.add_variable("lambda", FIRST, SCALAR);
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()
146  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
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 }
This is the EquationSystems class.
ConstFunction that simply returns 0.
Definition: zero_function.h:35
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.
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
MeshBase & mesh
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
SolverPackage default_solver_package()
Definition: libmesh.C:995
const DofMap & get_dof_map() const
Definition: system.h:2030
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
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
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 add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
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 assemble_elasticity(EquationSystems &es, const std::string &system_name)