libMesh
Public Member Functions | Private Attributes | List of all members
LinearElasticity Class Reference
Inheritance diagram for LinearElasticity:
[legend]

Public Member Functions

 LinearElasticity (EquationSystems &es_in)
 
Real kronecker_delta (unsigned int i, unsigned int j)
 Kronecker delta function. More...
 
Real elasticity_tensor (unsigned int i, unsigned int j, unsigned int k, unsigned int l)
 Evaluate the fourth order tensor (C_ijkl) that relates stress to strain. More...
 
void assemble ()
 Assemble the system matrix and right-hand side vector. More...
 
void compute_stresses ()
 

Private Attributes

EquationSystemses
 

Detailed Description

Definition at line 113 of file systems_of_equations_ex6.C.

Constructor & Destructor Documentation

LinearElasticity::LinearElasticity ( EquationSystems es_in)

Definition at line 120 of file systems_of_equations_ex6.C.

120  :
121  es(es_in)
122  {}

Member Function Documentation

void LinearElasticity::assemble ( )
virtual

Assemble the system matrix and right-hand side vector.

Implements libMesh::System::Assembly.

Definition at line 156 of file systems_of_equations_ex6.C.

References libMesh::MeshBase::active_local_element_ptr_range(), 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(), elasticity_tensor(), libMesh::MeshBase::get_boundary_info(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::BoundaryInfo::has_boundary_id(), libmesh_nullptr, libMesh::ImplicitSystem::matrix, mesh, libMesh::MeshBase::mesh_dimension(), libMesh::DenseMatrix< T >::resize(), libMesh::ExplicitSystem::rhs, side, libMesh::DenseSubVector< T >::size(), libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

157  {
158  const MeshBase & mesh = es.get_mesh();
159 
160  const unsigned int dim = mesh.mesh_dimension();
161 
162  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Elasticity");
163 
164  const unsigned int u_var = system.variable_number ("u");
165 
166  const DofMap & dof_map = system.get_dof_map();
167  FEType fe_type = dof_map.variable_type(u_var);
168  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
169  QGauss qrule (dim, fe_type.default_quadrature_order());
170  fe->attach_quadrature_rule (&qrule);
171 
172  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
173  QGauss qface(dim-1, fe_type.default_quadrature_order());
174  fe_face->attach_quadrature_rule (&qface);
175 
176  const std::vector<Real> & JxW = fe->get_JxW();
177  const std::vector<std::vector<Real>> & phi = fe->get_phi();
178  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
179 
181  DenseSubMatrix<Number> Ke_var[3][3] =
182  {
186  };
187 
189 
190  DenseSubVector<Number> Fe_var[3] =
194 
195  std::vector<dof_id_type> dof_indices;
196  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
197 
198  for (const auto & elem : mesh.active_local_element_ptr_range())
199  {
200  dof_map.dof_indices (elem, dof_indices);
201  for (unsigned int var=0; var<3; var++)
202  dof_map.dof_indices (elem, dof_indices_var[var], var);
203 
204  const unsigned int n_dofs = dof_indices.size();
205  const unsigned int n_var_dofs = dof_indices_var[0].size();
206 
207  fe->reinit (elem);
208 
209  Ke.resize (n_dofs, n_dofs);
210  for (unsigned int var_i=0; var_i<3; var_i++)
211  for (unsigned int var_j=0; var_j<3; var_j++)
212  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
213 
214  Fe.resize (n_dofs);
215  for (unsigned int var=0; var<3; var++)
216  Fe_var[var].reposition (var*n_var_dofs, n_var_dofs);
217 
218  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
219  {
220  // assemble \int_Omega C_ijkl u_k,l v_i,j \dx
221  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
222  for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
223  for (unsigned int i=0; i<3; i++)
224  for (unsigned int j=0; j<3; j++)
225  for (unsigned int k=0; k<3; k++)
226  for (unsigned int l=0; l<3; l++)
227  Ke_var[i][k](dof_i,dof_j) +=
228  JxW[qp] * elasticity_tensor(i,j,k,l) * dphi[dof_j][qp](l) * dphi[dof_i][qp](j);
229 
230  // assemble \int_Omega f_i v_i \dx
231  DenseVector<Number> f_vec(3);
232  f_vec(0) = 0.;
233  f_vec(1) = 0.;
234  f_vec(2) = -1.;
235  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
236  for (unsigned int i=0; i<3; i++)
237  Fe_var[i](dof_i) += JxW[qp] * (f_vec(i) * phi[dof_i][qp]);
238  }
239 
240  // assemble \int_\Gamma g_i v_i \ds
241  DenseVector<Number> g_vec(3);
242  g_vec(0) = 0.;
243  g_vec(1) = 0.;
244  g_vec(2) = -1.;
245  {
246  for (auto side : elem->side_index_range())
247  if (elem->neighbor_ptr(side) == libmesh_nullptr)
248  {
249  const std::vector<std::vector<Real>> & phi_face = fe_face->get_phi();
250  const std::vector<Real> & JxW_face = fe_face->get_JxW();
251 
252  fe_face->reinit(elem, side);
253 
254  // Apply a traction
255  for (unsigned int qp=0; qp<qface.n_points(); qp++)
256  if (mesh.get_boundary_info().has_boundary_id(elem, side, BOUNDARY_ID_MAX_X))
257  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
258  for (unsigned int i=0; i<3; i++)
259  Fe_var[i](dof_i) += JxW_face[qp] * (g_vec(i) * phi_face[dof_i][qp]);
260  }
261  }
262 
263  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
264 
265  system.matrix->add_matrix (Ke, dof_indices);
266  system.rhs->add_vector (Fe, dof_indices);
267  }
268  }
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.
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.
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.
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
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
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
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
const T_sys & get_system(const std::string &name) const
virtual unsigned int size() const libmesh_override
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 LinearElasticity::compute_stresses ( )

Definition at line 271 of file systems_of_equations_ex6.C.

References libMesh::MeshBase::active_local_element_ptr_range(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DenseMatrix< T >::cholesky_solve(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), elasticity_tensor(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::System::n_vars(), std::pow(), libMesh::System::solution, libMesh::System::update(), libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

Referenced by main().

272  {
273  const MeshBase & mesh = es.get_mesh();
274  const unsigned int dim = mesh.mesh_dimension();
275 
276  LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>("Elasticity");
277 
278  unsigned int displacement_vars[3];
279  displacement_vars[0] = system.variable_number ("u");
280  displacement_vars[1] = system.variable_number ("v");
281  displacement_vars[2] = system.variable_number ("w");
282  const unsigned int u_var = system.variable_number ("u");
283 
284  const DofMap & dof_map = system.get_dof_map();
285  FEType fe_type = dof_map.variable_type(u_var);
286  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
287  QGauss qrule (dim, fe_type.default_quadrature_order());
288  fe->attach_quadrature_rule (&qrule);
289 
290  const std::vector<Real> & JxW = fe->get_JxW();
291  const std::vector<std::vector<Real>> & phi = fe->get_phi();
292  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
293 
294  // Also, get a reference to the ExplicitSystem
295  ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
296  const DofMap & stress_dof_map = stress_system.get_dof_map();
297  unsigned int sigma_vars[6];
298  sigma_vars[0] = stress_system.variable_number ("sigma_00");
299  sigma_vars[1] = stress_system.variable_number ("sigma_01");
300  sigma_vars[2] = stress_system.variable_number ("sigma_02");
301  sigma_vars[3] = stress_system.variable_number ("sigma_11");
302  sigma_vars[4] = stress_system.variable_number ("sigma_12");
303  sigma_vars[5] = stress_system.variable_number ("sigma_22");
304  unsigned int vonMises_var = stress_system.variable_number ("vonMises");
305 
306  // Storage for the stress dof indices on each element
307  std::vector<std::vector<dof_id_type>> dof_indices_var(system.n_vars());
308  std::vector<dof_id_type> stress_dof_indices_var;
309  std::vector<dof_id_type> vonmises_dof_indices_var;
310 
311  for (const auto & elem : mesh.active_local_element_ptr_range())
312  {
313  for (unsigned int var=0; var<3; var++)
314  dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
315 
316  const unsigned int n_var_dofs = dof_indices_var[0].size();
317 
318  fe->reinit (elem);
319 
320  std::vector<DenseMatrix<Number>> stress_tensor_qp(qrule.n_points());
321  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
322  {
323  stress_tensor_qp[qp].resize(3,3);
324 
325  // Row is variable u1, u2, or u3, column is x, y, or z
326  DenseMatrix<Number> grad_u(3,3);
327  for (unsigned int var_i=0; var_i<3; var_i++)
328  for (unsigned int var_j=0; var_j<3; var_j++)
329  for (unsigned int j=0; j<n_var_dofs; j++)
330  grad_u(var_i,var_j) += dphi[j][qp](var_j) * system.current_solution(dof_indices_var[var_i][j]);
331 
332  for (unsigned int var_i=0; var_i<3; var_i++)
333  for (unsigned int var_j=0; var_j<3; var_j++)
334  for (unsigned int k=0; k<3; k++)
335  for (unsigned int l=0; l<3; l++)
336  stress_tensor_qp[qp](var_i,var_j) += elasticity_tensor(var_i,var_j,k,l) * grad_u(k,l);
337  }
338 
339  stress_dof_map.dof_indices (elem, vonmises_dof_indices_var, vonMises_var);
340  std::vector<DenseMatrix<Number>> elem_sigma_vec(vonmises_dof_indices_var.size());
341  for (std::size_t index=0; index<elem_sigma_vec.size(); index++)
342  elem_sigma_vec[index].resize(3,3);
343 
344  // Below we project each component of the stress tensor onto a L2_LAGRANGE discretization.
345  // Note that this gives a discontinuous stress plot on element boundaries, which is
346  // appropriate. We then also get the von Mises stress from the projected stress tensor.
347  unsigned int stress_var_index = 0;
348  for (unsigned int var_i=0; var_i<3; var_i++)
349  for (unsigned int var_j=var_i; var_j<3; var_j++)
350  {
351  stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[stress_var_index]);
352 
353  const unsigned int n_proj_dofs = stress_dof_indices_var.size();
354 
355  DenseMatrix<Real> Me(n_proj_dofs, n_proj_dofs);
356  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
357  {
358  for(unsigned int i=0; i<n_proj_dofs; i++)
359  for(unsigned int j=0; j<n_proj_dofs; j++)
360  {
361  Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
362  }
363  }
364 
365  DenseVector<Number> Fe(n_proj_dofs);
366  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
367  for(unsigned int i=0; i<n_proj_dofs; i++)
368  {
369  Fe(i) += JxW[qp] * stress_tensor_qp[qp](var_i,var_j) * phi[i][qp];
370  }
371 
372  DenseVector<Number> projected_data;
373  Me.cholesky_solve(Fe, projected_data);
374 
375  for(unsigned int index=0; index<n_proj_dofs; index++)
376  {
377  dof_id_type dof_index = stress_dof_indices_var[index];
378  if ((stress_system.solution->first_local_index() <= dof_index) &&
379  (dof_index < stress_system.solution->last_local_index()))
380  stress_system.solution->set(dof_index, projected_data(index));
381 
382  elem_sigma_vec[index](var_i,var_j) = projected_data(index);
383  }
384 
385  stress_var_index++;
386  }
387 
388  for (std::size_t index=0; index<elem_sigma_vec.size(); index++)
389  {
390  elem_sigma_vec[index](1,0) = elem_sigma_vec[index](0,1);
391  elem_sigma_vec[index](2,0) = elem_sigma_vec[index](0,2);
392  elem_sigma_vec[index](2,1) = elem_sigma_vec[index](1,2);
393 
394  // Get the von Mises stress from the projected stress tensor
395  Number vonMises_value = std::sqrt(0.5*(pow(elem_sigma_vec[index](0,0) - elem_sigma_vec[index](1,1), 2.) +
396  pow(elem_sigma_vec[index](1,1) - elem_sigma_vec[index](2,2), 2.) +
397  pow(elem_sigma_vec[index](2,2) - elem_sigma_vec[index](0,0), 2.) +
398  6.*(pow(elem_sigma_vec[index](0,1), 2.) +
399  pow(elem_sigma_vec[index](1,2), 2.) +
400  pow(elem_sigma_vec[index](2,0), 2.))));
401 
402  dof_id_type dof_index = vonmises_dof_indices_var[index];
403 
404  if ((stress_system.solution->first_local_index() <= dof_index) &&
405  (dof_index < stress_system.solution->last_local_index()))
406  stress_system.solution->set(dof_index, vonMises_value);
407  }
408  }
409 
410  // Should call close and update when we set vector entries directly
411  stress_system.solution->close();
412  stress_system.update();
413  }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
This class provides a specific system class.
MeshBase & mesh
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
const DofMap & get_dof_map() const
Definition: system.h:2030
Order default_quadrature_order() const
Definition: fe_type.h:332
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
double pow(double a, int b)
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
This class implements specific orders of Gauss quadrature.
const MeshBase & get_mesh() const
const T_sys & get_system(const std::string &name) const
unsigned int n_vars() const
Definition: system.h:2086
The ExplicitSystem provides only "right hand side" storage, which should be sufficient for solving mo...
uint8_t dof_id_type
Definition: id_types.h:64
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 LinearElasticity::elasticity_tensor ( unsigned int  i,
unsigned int  j,
unsigned int  k,
unsigned int  l 
)

Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.

Definition at line 136 of file systems_of_equations_ex6.C.

References kronecker_delta(), and libMesh::Real.

140  {
141  // Hard code material parameters for the sake of simplicity
142  const Real poisson_ratio = 0.3;
143  const Real young_modulus = 1.;
144 
145  // Define the Lame constants
146  const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
147  const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
148 
149  return lambda_1 * kronecker_delta(i, j) * kronecker_delta(k, l) +
150  lambda_2 * (kronecker_delta(i, k) * kronecker_delta(j, l) + kronecker_delta(i, l) * kronecker_delta(j, k));
151  }
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real kronecker_delta(unsigned int i, unsigned int j)
Kronecker delta function.
Real LinearElasticity::kronecker_delta ( unsigned int  i,
unsigned int  j 
)

Kronecker delta function.

Definition at line 127 of file systems_of_equations_ex6.C.

129  {
130  return i == j ? 1. : 0.;
131  }

Member Data Documentation

EquationSystems& LinearElasticity::es
private

Definition at line 116 of file systems_of_equations_ex6.C.


The documentation for this class was generated from the following file: