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

Public Member Functions

 LargeDeformationElasticity (EquationSystems &es_in)
 
Real kronecker_delta (unsigned int i, unsigned int j)
 Kronecker delta function. More...
 
Real elasticity_tensor (Real young_modulus, Real poisson_ratio, 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...
 
virtual void jacobian (const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, NonlinearImplicitSystem &)
 Evaluate the Jacobian of the nonlinear system. More...
 
virtual void residual (const NumericVector< Number > &soln, NumericVector< Number > &residual, NonlinearImplicitSystem &)
 Evaluate the residual of the nonlinear system. More...
 
void compute_stresses ()
 Compute the Cauchy stress for the current solution. More...
 
virtual void residual (const NumericVector< Number > &X, NumericVector< Number > &R, sys_type &S)=0
 Residual function. More...
 
virtual void jacobian (const NumericVector< Number > &X, SparseMatrix< Number > &J, sys_type &S)=0
 Jacobian function. More...
 

Private Attributes

EquationSystemses
 

Detailed Description

Definition at line 84 of file systems_of_equations_ex7.C.

Constructor & Destructor Documentation

LargeDeformationElasticity::LargeDeformationElasticity ( EquationSystems es_in)

Definition at line 92 of file systems_of_equations_ex7.C.

92  :
93  es(es_in)
94  {}

Member Function Documentation

void LargeDeformationElasticity::compute_stresses ( )

Compute the Cauchy stress for the current solution.

Definition at line 392 of file systems_of_equations_ex7.C.

References libMesh::DenseMatrix< T >::add(), libMesh::FEGenericBase< OutputType >::build(), libMesh::System::current_solution(), libMesh::FEType::default_quadrature_order(), libMesh::DenseMatrix< T >::det(), dim, libMesh::DofMap::dof_indices(), elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), libMesh::DenseMatrix< T >::left_multiply(), mesh, libMesh::System::n_vars(), libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseMatrix< T >::resize(), libMesh::DenseMatrix< T >::right_multiply_transpose(), libMesh::DenseMatrix< T >::scale(), libMesh::System::solution, libMesh::System::update(), libMesh::System::variable_number(), and libMesh::DofMap::variable_type().

Referenced by main().

393  {
394  const Real young_modulus = es.parameters.get<Real>("young_modulus");
395  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
396 
397  const MeshBase & mesh = es.get_mesh();
398  const unsigned int dim = mesh.mesh_dimension();
399 
400  NonlinearImplicitSystem & system =
401  es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");
402 
403  unsigned int displacement_vars[3];
404  displacement_vars[0] = system.variable_number ("u");
405  displacement_vars[1] = system.variable_number ("v");
406  displacement_vars[2] = system.variable_number ("w");
407  const unsigned int u_var = system.variable_number ("u");
408 
409  const DofMap & dof_map = system.get_dof_map();
410  FEType fe_type = dof_map.variable_type(u_var);
411  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
412  QGauss qrule (dim, fe_type.default_quadrature_order());
413  fe->attach_quadrature_rule (&qrule);
414 
415  const std::vector<Real> & JxW = fe->get_JxW();
416  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
417 
418  // Also, get a reference to the ExplicitSystem
419  ExplicitSystem & stress_system = es.get_system<ExplicitSystem>("StressSystem");
420  const DofMap & stress_dof_map = stress_system.get_dof_map();
421  unsigned int sigma_vars[6];
422  sigma_vars[0] = stress_system.variable_number ("sigma_00");
423  sigma_vars[1] = stress_system.variable_number ("sigma_01");
424  sigma_vars[2] = stress_system.variable_number ("sigma_02");
425  sigma_vars[3] = stress_system.variable_number ("sigma_11");
426  sigma_vars[4] = stress_system.variable_number ("sigma_12");
427  sigma_vars[5] = stress_system.variable_number ("sigma_22");
428 
429  // Storage for the stress dof indices on each element
430  std::vector<std::vector<dof_id_type>> dof_indices_var(system.n_vars());
431  std::vector<dof_id_type> stress_dof_indices_var;
432 
433  // To store the stress tensor on each element
434  DenseMatrix<Number> elem_avg_stress_tensor(3, 3);
435 
436  for (const auto & elem : mesh.active_local_element_ptr_range())
437  {
438  for (unsigned int var=0; var<3; var++)
439  dof_map.dof_indices (elem, dof_indices_var[var], displacement_vars[var]);
440 
441  const unsigned int n_var_dofs = dof_indices_var[0].size();
442 
443  fe->reinit (elem);
444 
445  // clear the stress tensor
446  elem_avg_stress_tensor.resize(3, 3);
447 
448  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
449  {
450  DenseMatrix<Number> grad_u(3, 3);
451  // Row is variable u1, u2, or u3, column is x, y, or z
452  for (unsigned int var_i=0; var_i<3; var_i++)
453  for (unsigned int var_j=0; var_j<3; var_j++)
454  for (unsigned int j=0; j<n_var_dofs; j++)
455  grad_u(var_i,var_j) += dphi[j][qp](var_j) * system.current_solution(dof_indices_var[var_i][j]);
456 
457  DenseMatrix<Number> strain_tensor(3, 3);
458  for (unsigned int i=0; i<3; i++)
459  for (unsigned int j=0; j<3; j++)
460  {
461  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));
462 
463  for (unsigned int k=0; k<3; k++)
464  strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
465  }
466 
467  // Define the deformation gradient
468  DenseMatrix<Number> F(3, 3);
469  F = grad_u;
470  for (unsigned int var=0; var<3; var++)
471  F(var, var) += 1.;
472 
473  DenseMatrix<Number> stress_tensor(3, 3);
474  for (unsigned int i=0; i<3; i++)
475  for (unsigned int j=0; j<3; j++)
476  for (unsigned int k=0; k<3; k++)
477  for (unsigned int l=0; l<3; l++)
478  stress_tensor(i,j) +=
479  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k, l);
480 
481  // stress_tensor now holds the second Piola-Kirchoff stress (PK2) at point qp.
482  // However, in this example we want to compute the Cauchy stress which is given by
483  // 1/det(F) * F * PK2 * F^T, hence we now apply this transformation.
484  stress_tensor.scale(1./F.det());
485  stress_tensor.left_multiply(F);
486  stress_tensor.right_multiply_transpose(F);
487 
488  // We want to plot the average Cauchy stress on each element, hence
489  // we integrate stress_tensor
490  elem_avg_stress_tensor.add(JxW[qp], stress_tensor);
491  }
492 
493  // Get the average stress per element by dividing by volume
494  elem_avg_stress_tensor.scale(1./elem->volume());
495 
496  // load elem_sigma data into stress_system
497  unsigned int stress_var_index = 0;
498  for (unsigned int i=0; i<3; i++)
499  for (unsigned int j=i; j<3; j++)
500  {
501  stress_dof_map.dof_indices (elem, stress_dof_indices_var, sigma_vars[stress_var_index]);
502 
503  // We are using CONSTANT MONOMIAL basis functions, hence we only need to get
504  // one dof index per variable
505  dof_id_type dof_index = stress_dof_indices_var[0];
506 
507  if ((stress_system.solution->first_local_index() <= dof_index) &&
508  (dof_index < stress_system.solution->last_local_index()))
509  stress_system.solution->set(dof_index, elem_avg_stress_tensor(i,j));
510 
511  stress_var_index++;
512  }
513  }
514 
515  // Should call close and update when we set vector entries directly
516  stress_system.solution->close();
517  stress_system.update();
518  }
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
MeshBase & mesh
const T & get(const std::string &) const
Definition: parameters.h:430
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
This class provides a specific system class.
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:192
Real elasticity_tensor(Real young_modulus, Real poisson_ratio, unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
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 LargeDeformationElasticity::elasticity_tensor ( Real  young_modulus,
Real  poisson_ratio,
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 108 of file systems_of_equations_ex7.C.

References kronecker_delta(), and libMesh::Real.

114  {
115  // Define the Lame constants
116  const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
117  const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
118 
119  return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l) +
120  lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l) + kronecker_delta(i,l) * kronecker_delta(j,k));
121  }
Real kronecker_delta(unsigned int i, unsigned int j)
Kronecker delta function.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void libMesh::NonlinearImplicitSystem::ComputeJacobian::jacobian ( const NumericVector< Number > &  X,
SparseMatrix< Number > &  J,
sys_type S 
)
pure virtualinherited

Jacobian function.

This function will be called to compute the jacobian and must be implemented by the user in a derived class.

Referenced by libMesh::__libmesh_petsc_snes_jacobian(), libMesh::Problem_Interface::computeJacobian(), and libMesh::Problem_Interface::computePreconditioner().

virtual void LargeDeformationElasticity::jacobian ( const NumericVector< Number > &  soln,
SparseMatrix< Number > &  jacobian,
NonlinearImplicitSystem  
)
virtual

Evaluate the Jacobian of the nonlinear system.

Definition at line 127 of file systems_of_equations_ex7.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DofMap::constrain_element_matrix(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseMatrix< T >::resize(), libMesh::System::variable_number(), libMesh::DofMap::variable_type(), and libMesh::SparseMatrix< T >::zero().

130  {
131  const Real young_modulus = es.parameters.get<Real>("young_modulus");
132  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
133 
134  const MeshBase & mesh = es.get_mesh();
135  const unsigned int dim = mesh.mesh_dimension();
136 
137  NonlinearImplicitSystem & system =
138  es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");
139 
140  const unsigned int u_var = system.variable_number ("u");
141 
142  const DofMap & dof_map = system.get_dof_map();
143 
144  FEType fe_type = dof_map.variable_type(u_var);
145  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
146  QGauss qrule (dim, fe_type.default_quadrature_order());
147  fe->attach_quadrature_rule (&qrule);
148 
149  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
150  QGauss qface (dim-1, fe_type.default_quadrature_order());
151  fe_face->attach_quadrature_rule (&qface);
152 
153  const std::vector<Real> & JxW = fe->get_JxW();
154  const std::vector<std::vector<Real>> & phi = fe->get_phi();
155  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
156 
158  DenseSubMatrix<Number> Ke_var[3][3] =
159  {
163  };
164 
165  std::vector<dof_id_type> dof_indices;
166  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
167 
168  jacobian.zero();
169 
170  for (const auto & elem : mesh.active_local_element_ptr_range())
171  {
172  dof_map.dof_indices (elem, dof_indices);
173  for (unsigned int var=0; var<3; var++)
174  dof_map.dof_indices (elem, dof_indices_var[var], var);
175 
176  const unsigned int n_dofs = dof_indices.size();
177  const unsigned int n_var_dofs = dof_indices_var[0].size();
178 
179  fe->reinit (elem);
180 
181  Ke.resize (n_dofs, n_dofs);
182  for (unsigned int var_i=0; var_i<3; var_i++)
183  for (unsigned int var_j=0; var_j<3; var_j++)
184  Ke_var[var_i][var_j].reposition (var_i*n_var_dofs, var_j*n_var_dofs, n_var_dofs, n_var_dofs);
185 
186  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
187  {
188  DenseVector<Number> u_vec(3);
189  DenseMatrix<Number> grad_u(3, 3);
190  for (unsigned int var_i=0; var_i<3; var_i++)
191  {
192  for (unsigned int j=0; j<n_var_dofs; j++)
193  u_vec(var_i) += phi[j][qp]*soln(dof_indices_var[var_i][j]);
194 
195  // Row is variable u1, u2, or u3, column is x, y, or z
196  for (unsigned int var_j=0; var_j<3; var_j++)
197  for (unsigned int j=0; j<n_var_dofs; j++)
198  grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
199  }
200 
201  DenseMatrix<Number> strain_tensor(3, 3);
202  for (unsigned int i=0; i<3; i++)
203  for (unsigned int j=0; j<3; j++)
204  {
205  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));
206 
207  for (unsigned int k=0; k<3; k++)
208  strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
209  }
210 
211  // Define the deformation gradient
212  DenseMatrix<Number> F(3, 3);
213  F = grad_u;
214  for (unsigned int var=0; var<3; var++)
215  F(var, var) += 1.;
216 
217  DenseMatrix<Number> stress_tensor(3, 3);
218 
219  for (unsigned int i=0; i<3; i++)
220  for (unsigned int j=0; j<3; j++)
221  for (unsigned int k=0; k<3; k++)
222  for (unsigned int l=0; l<3; l++)
223  stress_tensor(i,j) +=
224  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k, l);
225 
226  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
227  for (unsigned int dof_j=0; dof_j<n_var_dofs; dof_j++)
228  {
229  for (unsigned int i=0; i<3; i++)
230  for (unsigned int j=0; j<3; j++)
231  for (unsigned int m=0; m<3; m++)
232  Ke_var[i][i](dof_i,dof_j) += JxW[qp] *
233  (-dphi[dof_j][qp](m) * stress_tensor(m,j) * dphi[dof_i][qp](j));
234 
235  for (unsigned int i=0; i<3; i++)
236  for (unsigned int j=0; j<3; j++)
237  for (unsigned int k=0; k<3; k++)
238  for (unsigned int l=0; l<3; l++)
239  {
240  Number FxC_ijkl = 0.;
241  for (unsigned int m=0; m<3; m++)
242  FxC_ijkl += F(i,m) * elasticity_tensor(young_modulus, poisson_ratio, m, j, k, l);
243 
244  Ke_var[i][k](dof_i,dof_j) += JxW[qp] *
245  (-0.5 * FxC_ijkl * dphi[dof_j][qp](l) * dphi[dof_i][qp](j));
246 
247  Ke_var[i][l](dof_i,dof_j) += JxW[qp] *
248  (-0.5 * FxC_ijkl * dphi[dof_j][qp](k) * dphi[dof_i][qp](j));
249 
250  for (unsigned int n=0; n<3; n++)
251  Ke_var[i][n](dof_i,dof_j) += JxW[qp] *
252  (-0.5 * FxC_ijkl * (dphi[dof_j][qp](k) * grad_u(n,l) + dphi[dof_j][qp](l) * grad_u(n,k)) * dphi[dof_i][qp](j));
253  }
254  }
255  }
256 
257  dof_map.constrain_element_matrix (Ke, dof_indices);
258  jacobian.add_matrix (Ke, dof_indices);
259  }
260  }
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
MeshBase & mesh
const T & get(const std::string &) const
Definition: parameters.h:430
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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
virtual void zero()=0
Set all entries to 0.
Order default_quadrature_order() const
Definition: fe_type.h:332
Defines a dense submatrix for use in Finite Element-type computations.
This class provides a specific system class.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void constrain_element_matrix(DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix.
Definition: dof_map.h:1793
Real elasticity_tensor(Real young_modulus, Real poisson_ratio, unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
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.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
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 LargeDeformationElasticity::kronecker_delta ( unsigned int  i,
unsigned int  j 
)

Kronecker delta function.

Definition at line 99 of file systems_of_equations_ex7.C.

101  {
102  return i == j ? 1. : 0.;
103  }
virtual void libMesh::NonlinearImplicitSystem::ComputeResidual::residual ( const NumericVector< Number > &  X,
NumericVector< Number > &  R,
sys_type S 
)
pure virtualinherited

Residual function.

This function will be called to compute the residual and must be implemented by the user in a derived class.

Referenced by libMesh::__libmesh_petsc_snes_residual(), libMesh::Problem_Interface::computeF(), and libMesh::NonlinearImplicitSystem::ComputeResidual::~ComputeResidual().

virtual void LargeDeformationElasticity::residual ( const NumericVector< Number > &  soln,
NumericVector< Number > &  residual,
NonlinearImplicitSystem  
)
virtual

Evaluate the residual of the nonlinear system.

Definition at line 265 of file systems_of_equations_ex7.C.

References libMesh::NumericVector< T >::add_vector(), libMesh::FEGenericBase< OutputType >::build(), libMesh::DofMap::constrain_element_vector(), libMesh::FEType::default_quadrature_order(), dim, libMesh::DofMap::dof_indices(), elasticity_tensor(), libMesh::Parameters::get(), libMesh::System::get_dof_map(), libMesh::EquationSystems::get_mesh(), libMesh::EquationSystems::get_system(), mesh, libMesh::EquationSystems::parameters, libMesh::Real, libMesh::DenseVector< T >::resize(), libMesh::DenseSubVector< T >::size(), libMesh::System::variable_number(), libMesh::DofMap::variable_type(), and libMesh::NumericVector< T >::zero().

268  {
269  const Real young_modulus = es.parameters.get<Real>("young_modulus");
270  const Real poisson_ratio = es.parameters.get<Real>("poisson_ratio");
271  const Real forcing_magnitude = es.parameters.get<Real>("forcing_magnitude");
272 
273  const MeshBase & mesh = es.get_mesh();
274  const unsigned int dim = mesh.mesh_dimension();
275 
276  NonlinearImplicitSystem & system =
277  es.get_system<NonlinearImplicitSystem>("NonlinearElasticity");
278 
279  const unsigned int u_var = system.variable_number ("u");
280 
281  const DofMap & dof_map = system.get_dof_map();
282 
283  FEType fe_type = dof_map.variable_type(u_var);
284  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
285  QGauss qrule (dim, fe_type.default_quadrature_order());
286  fe->attach_quadrature_rule (&qrule);
287 
288  UniquePtr<FEBase> fe_face (FEBase::build(dim, fe_type));
289  QGauss qface (dim-1, fe_type.default_quadrature_order());
290  fe_face->attach_quadrature_rule (&qface);
291 
292  const std::vector<Real> & JxW = fe->get_JxW();
293  const std::vector<std::vector<Real>> & phi = fe->get_phi();
294  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
295 
297 
298  DenseSubVector<Number> Re_var[3] =
302 
303  std::vector<dof_id_type> dof_indices;
304  std::vector<std::vector<dof_id_type>> dof_indices_var(3);
305 
306  residual.zero();
307 
308  for (const auto & elem : mesh.active_local_element_ptr_range())
309  {
310  dof_map.dof_indices (elem, dof_indices);
311  for (unsigned int var=0; var<3; var++)
312  dof_map.dof_indices (elem, dof_indices_var[var], var);
313 
314  const unsigned int n_dofs = dof_indices.size();
315  const unsigned int n_var_dofs = dof_indices_var[0].size();
316 
317  fe->reinit (elem);
318 
319  Re.resize (n_dofs);
320  for (unsigned int var=0; var<3; var++)
321  Re_var[var].reposition (var*n_var_dofs, n_var_dofs);
322 
323  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
324  {
325  DenseVector<Number> u_vec(3);
326  DenseMatrix<Number> grad_u(3, 3);
327  for (unsigned int var_i=0; var_i<3; var_i++)
328  {
329  for (unsigned int j=0; j<n_var_dofs; j++)
330  u_vec(var_i) += phi[j][qp]*soln(dof_indices_var[var_i][j]);
331 
332  // Row is variable u, v, or w column is x, y, or z
333  for (unsigned int var_j=0; var_j<3; var_j++)
334  for (unsigned int j=0; j<n_var_dofs; j++)
335  grad_u(var_i,var_j) += dphi[j][qp](var_j)*soln(dof_indices_var[var_i][j]);
336  }
337 
338  DenseMatrix<Number> strain_tensor(3, 3);
339  for (unsigned int i=0; i<3; i++)
340  for (unsigned int j=0; j<3; j++)
341  {
342  strain_tensor(i,j) += 0.5 * (grad_u(i,j) + grad_u(j,i));
343 
344  for (unsigned int k=0; k<3; k++)
345  strain_tensor(i,j) += 0.5 * grad_u(k,i)*grad_u(k,j);
346  }
347 
348  // Define the deformation gradient
349  DenseMatrix<Number> F(3, 3);
350  F = grad_u;
351  for (unsigned int var=0; var<3; var++)
352  F(var, var) += 1.;
353 
354  DenseMatrix<Number> stress_tensor(3, 3);
355 
356  for (unsigned int i=0; i<3; i++)
357  for (unsigned int j=0; j<3; j++)
358  for (unsigned int k=0; k<3; k++)
359  for (unsigned int l=0; l<3; l++)
360  stress_tensor(i,j) +=
361  elasticity_tensor(young_modulus, poisson_ratio, i, j, k, l) * strain_tensor(k,l);
362 
363  DenseVector<Number> f_vec(3);
364  f_vec(0) = 0.;
365  f_vec(1) = 0.;
366  f_vec(2) = -forcing_magnitude;
367 
368  for (unsigned int dof_i=0; dof_i<n_var_dofs; dof_i++)
369  for (unsigned int i=0; i<3; i++)
370  {
371  for (unsigned int j=0; j<3; j++)
372  {
373  Number FxStress_ij = 0.;
374  for (unsigned int m=0; m<3; m++)
375  FxStress_ij += F(i,m) * stress_tensor(m,j);
376 
377  Re_var[i](dof_i) += JxW[qp] * (-FxStress_ij * dphi[dof_i][qp](j));
378  }
379 
380  Re_var[i](dof_i) += JxW[qp] * (f_vec(i) * phi[dof_i][qp]);
381  }
382  }
383 
384  dof_map.constrain_element_vector (Re, dof_indices);
385  residual.add_vector (Re, dof_indices);
386  }
387  }
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
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]...
void constrain_element_vector(DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
Constrains the element vector.
Definition: dof_map.h:1802
MeshBase & mesh
const T & get(const std::string &) const
Definition: parameters.h:430
This is the MeshBase class.
Definition: mesh_base.h:68
virtual void zero()=0
Set all entries to zero.
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
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: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
This class provides a specific system class.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real elasticity_tensor(Real young_modulus, Real poisson_ratio, unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Evaluate the fourth order tensor (C_ijkl) that relates stress to strain.
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
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

Member Data Documentation

EquationSystems& LargeDeformationElasticity::es
private

Definition at line 88 of file systems_of_equations_ex7.C.


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