libMesh
Functions
systems_of_equations_ex4.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

◆ assemble_elasticity() [1/2]

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

Referenced by main().

◆ assemble_elasticity() [2/2]

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

Definition at line 164 of file systems_of_equations_ex4.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), 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::ImplicitSystem::get_system_matrix(), libMesh::BoundaryInfo::has_boundary_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::DenseSubVector< T >::reposition(), libMesh::DenseSubMatrix< T >::reposition(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

166 {
167  libmesh_assert_equal_to (system_name, "Elasticity");
168 
169  const MeshBase & mesh = es.get_mesh();
170 
171  const unsigned int dim = mesh.mesh_dimension();
172 
173  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Elasticity");
174 
175  const unsigned int u_var = system.variable_number ("u");
176  const unsigned int v_var = system.variable_number ("v");
177 
178  const DofMap & dof_map = system.get_dof_map();
179  FEType fe_type = dof_map.variable_type(0);
180  std::unique_ptr<FEBase> fe (FEBase::build(dim, fe_type));
181  QGauss qrule (dim, fe_type.default_quadrature_order());
182  fe->attach_quadrature_rule (&qrule);
183 
184  std::unique_ptr<FEBase> fe_face (FEBase::build(dim, fe_type));
185  QGauss qface(dim-1, fe_type.default_quadrature_order());
186  fe_face->attach_quadrature_rule (&qface);
187 
188  const std::vector<Real> & JxW = fe->get_JxW();
189  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
190 
193 
195  Kuu(Ke), Kuv(Ke),
196  Kvu(Ke), Kvv(Ke);
197 
199  Fu(Fe),
200  Fv(Fe);
201 
202  std::vector<dof_id_type> dof_indices;
203  std::vector<dof_id_type> dof_indices_u;
204  std::vector<dof_id_type> dof_indices_v;
205 
206  SparseMatrix<Number> & matrix = system.get_system_matrix();
207 
208  for (const auto & elem : mesh.active_local_element_ptr_range())
209  {
210  dof_map.dof_indices (elem, dof_indices);
211  dof_map.dof_indices (elem, dof_indices_u, u_var);
212  dof_map.dof_indices (elem, dof_indices_v, v_var);
213 
214  const unsigned int n_dofs = dof_indices.size();
215  const unsigned int n_u_dofs = dof_indices_u.size();
216  const unsigned int n_v_dofs = dof_indices_v.size();
217 
218  fe->reinit (elem);
219 
220  Ke.resize (n_dofs, n_dofs);
221  Fe.resize (n_dofs);
222 
223  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
224  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
225 
226  Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
227  Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
228 
229  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
230  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
231 
232  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
233  {
234  for (unsigned int i=0; i<n_u_dofs; i++)
235  for (unsigned int j=0; j<n_u_dofs; j++)
236  {
237  // Tensor indices
238  unsigned int C_i, C_j, C_k, C_l;
239  C_i=0, C_k=0;
240 
241  C_j=0, C_l=0;
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  C_j=1, C_l=0;
245  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));
246 
247  C_j=0, C_l=1;
248  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));
249 
250  C_j=1, C_l=1;
251  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));
252  }
253 
254  for (unsigned int i=0; i<n_u_dofs; i++)
255  for (unsigned int j=0; j<n_v_dofs; j++)
256  {
257  // Tensor indices
258  unsigned int C_i, C_j, C_k, C_l;
259  C_i=0, C_k=1;
260 
261  C_j=0, C_l=0;
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  C_j=1, C_l=0;
265  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));
266 
267  C_j=0, C_l=1;
268  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));
269 
270  C_j=1, C_l=1;
271  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));
272  }
273 
274  for (unsigned int i=0; i<n_v_dofs; i++)
275  for (unsigned int j=0; j<n_u_dofs; j++)
276  {
277  // Tensor indices
278  unsigned int C_i, C_j, C_k, C_l;
279  C_i=1, C_k=0;
280 
281  C_j=0, C_l=0;
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  C_j=1, C_l=0;
285  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));
286 
287  C_j=0, C_l=1;
288  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));
289 
290  C_j=1, C_l=1;
291  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));
292  }
293 
294  for (unsigned int i=0; i<n_v_dofs; i++)
295  for (unsigned int j=0; j<n_v_dofs; j++)
296  {
297  // Tensor indices
298  unsigned int C_i, C_j, C_k, C_l;
299  C_i=1, C_k=1;
300 
301  C_j=0, C_l=0;
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  C_j=1, C_l=0;
305  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));
306 
307  C_j=0, C_l=1;
308  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));
309 
310  C_j=1, C_l=1;
311  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));
312  }
313  }
314 
315  {
316  for (auto side : elem->side_index_range())
317  if (elem->neighbor_ptr(side) == nullptr)
318  {
319  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
320  const std::vector<Real> & JxW_face = fe_face->get_JxW();
321 
322  fe_face->reinit(elem, side);
323 
324  if (mesh.get_boundary_info().has_boundary_id (elem, side, 1)) // Apply a traction on the right side
325  {
326  for (unsigned int qp=0; qp<qface.n_points(); qp++)
327  for (unsigned int i=0; i<n_v_dofs; i++)
328  Fv(i) += JxW_face[qp] * (-1.) * phi_face[i][qp];
329  }
330  }
331  }
332 
333  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
334 
335  matrix.add_matrix (Ke, dof_indices);
336  system.rhs->add_vector (Fe, dof_indices);
337  }
338 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
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:2254
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:1992
unsigned int dim
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
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]...
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2144
MeshBase & mesh
NumericVector< Number > * rhs
The system matrix.
Order default_quadrature_order() const
Definition: fe_type.h:357
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:159
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
Definition: mesh_base.h:74
unsigned int variable_number(std::string_view var) const
Definition: system.C:1557
Real eval_elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Defines a dense subvector for use in finite element computations.
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:169
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.
Defines a dense submatrix for use in Finite Element-type computations.
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Definition: mesh_base.C:324
const DofMap & get_dof_map() const
Definition: system.h:2293
const SparseMatrix< Number > & get_system_matrix() const

◆ eval_elasticity_tensor()

Real eval_elasticity_tensor ( unsigned int  i,
unsigned int  j,
unsigned int  k,
unsigned int  l 
)

Definition at line 340 of file systems_of_equations_ex4.C.

References libMesh::Real.

Referenced by assemble_elasticity().

344 {
345  // Define the Poisson ratio
346  const Real nu = 0.3;
347 
348  // Define the Lame constants (lambda_1 and lambda_2) based on Poisson ratio
349  const Real lambda_1 = nu / ((1. + nu) * (1. - 2.*nu));
350  const Real lambda_2 = 0.5 / (1 + nu);
351 
352  // Define the Kronecker delta functions that we need here
353  Real delta_ij = (i == j) ? 1. : 0.;
354  Real delta_il = (i == l) ? 1. : 0.;
355  Real delta_ik = (i == k) ? 1. : 0.;
356  Real delta_jl = (j == l) ? 1. : 0.;
357  Real delta_jk = (j == k) ? 1. : 0.;
358  Real delta_kl = (k == l) ? 1. : 0.;
359 
360  return lambda_1 * delta_ij * delta_kl + lambda_2 * (delta_ik * delta_jl + delta_il * delta_jk);
361 }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 76 of file systems_of_equations_ex4.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::command_line_next(), libMesh::default_solver_package(), dim, 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::SECOND, libMesh::LinearImplicitSystem::solve(), and libMesh::MeshOutput< MT >::write_equation_systems().

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 
99  // Get the mesh size from the command line.
100  const int nx = libMesh::command_line_next("-nx", 50),
101  ny = libMesh::command_line_next("-ny", 10);
102 
104  nx, ny,
105  0., 1.,
106  0., 0.2,
107  QUAD9);
108 
109 
110  // Print information about the mesh to the screen.
111  mesh.print_info();
112 
113  // Create an equation systems object.
114  EquationSystems equation_systems (mesh);
115 
116  // Declare the system and its variables.
117  // Create a system named "Elasticity"
118  LinearImplicitSystem & system =
119  equation_systems.add_system<LinearImplicitSystem> ("Elasticity");
120 
122 
123 #ifdef LIBMESH_ENABLE_DIRICHLET
124  // Add two displacement variables, u and v, to the system
125  unsigned int u_var = system.add_variable("u", SECOND, LAGRANGE);
126  unsigned int v_var = system.add_variable("v", SECOND, LAGRANGE);
127 
128  // Create a ZeroFunction to initialize dirichlet_bc
129  ZeroFunction<> zf;
130 
131  // Construct a Dirichlet boundary condition object
132  // We impose a "clamped" boundary condition on the
133  // "left" boundary, i.e. bc_id = 3
134 
135  // Most DirichletBoundary users will want to supply a "locally
136  // indexed" functor
137  DirichletBoundary dirichlet_bc({3}, {u_var, v_var}, zf,
139 
140  // We must add the Dirichlet boundary condition _before_
141  // we call equation_systems.init()
142  system.get_dof_map().add_dirichlet_boundary(dirichlet_bc);
143 #endif // LIBMESH_ENABLE_DIRICHLET
144 
145  // Initialize the data structures for the equation system.
146  equation_systems.init();
147 
148  // Print information about the system to the screen.
149  equation_systems.print_info();
150 
151  // Solve the system
152  system.solve();
153 
154  // Plot the solution
155 #ifdef LIBMESH_HAVE_EXODUS_API
156  ExodusII_IO (mesh).write_equation_systems("displacement.e", equation_systems);
157 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
158 
159  // All done.
160  return 0;
161 }
T command_line_next(std::string name, T default_value)
Use GetPot&#39;s search()/next() functions to get following arguments from the command line...
Definition: libmesh.C:1011
This is the EquationSystems class.
ConstFunction that simply returns 0.
Definition: zero_function.h:38
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=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
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
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:90
virtual void solve() override
Assembles & solves the linear system A*x=b.
SolverPackage default_solver_package()
Definition: libmesh.C:1050
void assemble_elasticity(EquationSystems &es, const std::string &system_name)
void print_info(std::ostream &os=libMesh::out, const unsigned int verbosity=0, const bool global=true) const
Prints relevant information about the mesh.
Definition: mesh_base.C:1489
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1305
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:2109
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
const DofMap & get_dof_map() const
Definition: system.h:2293