libMesh
elasticity_system.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
6 // License as published by the Free Software Foundation; either
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 #include "elasticity_system.h"
19 
20 #include "libmesh/dof_map.h"
21 #include "libmesh/fe_base.h"
22 #include "libmesh/fem_context.h"
23 #include "libmesh/zero_function.h"
24 #include "libmesh/dirichlet_boundaries.h"
25 #include "libmesh/quadrature.h"
26 #include "libmesh/unsteady_solver.h"
27 
28 using namespace libMesh;
29 
31 {
32  _u_var = this->add_variable ("Ux", FIRST, LAGRANGE);
33  _v_var = this->add_variable ("Uy", FIRST, LAGRANGE);
34  _w_var = this->add_variable ("Uz", FIRST, LAGRANGE);
35 
36  this->time_evolving(_u_var,2);
37  this->time_evolving(_v_var,2);
38  this->time_evolving(_w_var,2);
39 
40  std::set<boundary_id_type> boundary_ids;
41  boundary_ids.insert(BOUNDARY_ID_MIN_X);
42  boundary_ids.insert(NODE_BOUNDARY_ID);
43  boundary_ids.insert(EDGE_BOUNDARY_ID);
44 
45  std::vector<unsigned int> variables;
46  variables.push_back(_u_var);
47  variables.push_back(_v_var);
48  variables.push_back(_w_var);
49 
50  ZeroFunction<> zf;
51 
52  // Most DirichletBoundary users will want to supply a "locally
53  // indexed" functor
54  DirichletBoundary dirichlet_bc(boundary_ids, variables, zf,
56 
57  this->get_dof_map().add_dirichlet_boundary(dirichlet_bc);
58 
59  // Do the parent's initialization after variables and boundary constraints are defined
61 }
62 
64 {
65  FEMContext & c = cast_ref<FEMContext &>(context);
66 
67  FEBase * u_elem_fe;
68  FEBase * u_side_fe;
69 
70  c.get_element_fe(_u_var, u_elem_fe);
71  c.get_side_fe(_u_var, u_side_fe);
72 
73  // We should prerequest all the data
74  // we will need to build the residuals.
75  u_elem_fe->get_JxW();
76  u_elem_fe->get_phi();
77  u_elem_fe->get_dphi();
78 
79  u_side_fe->get_JxW();
80  u_side_fe->get_phi();
81 }
82 
83 bool ElasticitySystem::element_time_derivative(bool request_jacobian,
84  DiffContext & context)
85 {
86  FEMContext & c = cast_ref<FEMContext &>(context);
87 
88  // If we have an unsteady solver, then we need to extract the corresponding
89  // velocity variable. This allows us to use either a FirstOrderUnsteadySolver
90  // or a SecondOrderUnsteadySolver. That is, we get back the velocity variable
91  // index for FirstOrderUnsteadySolvers or, if it's a SecondOrderUnsteadySolver,
92  // this is actually just giving us back the same variable index.
93 
94  // If we only wanted to use a SecondOrderUnsteadySolver, then this
95  // step would be unnecessary and we would just
96  // populate the _u_var, etc. blocks of the residual and Jacobian.
97  unsigned int u_dot_var = this->get_second_order_dot_var(_u_var);
98  unsigned int v_dot_var = this->get_second_order_dot_var(_v_var);
99  unsigned int w_dot_var = this->get_second_order_dot_var(_w_var);
100 
101  FEBase * u_elem_fe;
102  c.get_element_fe(_u_var, u_elem_fe);
103 
104  // The number of local degrees of freedom in each variable
105  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
106 
107  // Element Jacobian * quadrature weights for interior integration
108  const std::vector<Real> & JxW = u_elem_fe->get_JxW();
109 
110  const std::vector<std::vector<Real>> & phi = u_elem_fe->get_phi();
111  const std::vector<std::vector<RealGradient>> & grad_phi = u_elem_fe->get_dphi();
112 
113  DenseSubVector<Number> & Fu = c.get_elem_residual(u_dot_var);
114  DenseSubVector<Number> & Fv = c.get_elem_residual(v_dot_var);
115  DenseSubVector<Number> & Fw = c.get_elem_residual(w_dot_var);
116 
117  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_dot_var, _u_var);
118  DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_dot_var, _v_var);
119  DenseSubMatrix<Number> & Kww = c.get_elem_jacobian(w_dot_var, _w_var);
120  DenseSubMatrix<Number> & Kuv = c.get_elem_jacobian(u_dot_var, _v_var);
121  DenseSubMatrix<Number> & Kuw = c.get_elem_jacobian(u_dot_var, _w_var);
122  DenseSubMatrix<Number> & Kvu = c.get_elem_jacobian(v_dot_var, _u_var);
123  DenseSubMatrix<Number> & Kvw = c.get_elem_jacobian(v_dot_var, _w_var);
124  DenseSubMatrix<Number> & Kwu = c.get_elem_jacobian(w_dot_var, _u_var);
125  DenseSubMatrix<Number> & Kwv = c.get_elem_jacobian(w_dot_var, _v_var);
126 
127  unsigned int n_qpoints = c.get_element_qrule().n_points();
128 
129  Gradient body_force(0.0, 0.0, -1.0);
130 
131  for (unsigned int qp=0; qp != n_qpoints; qp++)
132  {
133  Gradient grad_u, grad_v, grad_w;
134  c.interior_gradient(_u_var, qp, grad_u);
135  c.interior_gradient(_v_var, qp, grad_v);
136  c.interior_gradient(_w_var, qp, grad_w);
137 
138  // Convenience
139  Tensor grad_U (grad_u, grad_v, grad_w);
140 
141  Tensor tau;
142  for (unsigned int i = 0; i < 3; i++)
143  for (unsigned int j = 0; j < 3; j++)
144  for (unsigned int k = 0; k < 3; k++)
145  for (unsigned int l = 0; l < 3; l++)
146  tau(i,j) += elasticity_tensor(i,j,k,l)*grad_U(k,l);
147 
148  for (unsigned int i=0; i != n_u_dofs; i++)
149  {
150  for (unsigned int alpha = 0; alpha < 3; alpha++)
151  {
152  Fu(i) += (tau(0,alpha)*grad_phi[i][qp](alpha) - body_force(0)*phi[i][qp])*JxW[qp];
153  Fv(i) += (tau(1,alpha)*grad_phi[i][qp](alpha) - body_force(1)*phi[i][qp])*JxW[qp];
154  Fw(i) += (tau(2,alpha)*grad_phi[i][qp](alpha) - body_force(2)*phi[i][qp])*JxW[qp];
155 
156  if (request_jacobian)
157  {
158  for (unsigned int j=0; j != n_u_dofs; j++)
159  {
160  for (unsigned int beta = 0; beta < 3; beta++)
161  {
162  // Convenience
163  const Real c0 = grad_phi[j][qp](beta)*c.get_elem_solution_derivative();
164 
165  Real dtau_uu = elasticity_tensor(0, alpha, 0, beta)*c0;
166  Real dtau_uv = elasticity_tensor(0, alpha, 1, beta)*c0;
167  Real dtau_uw = elasticity_tensor(0, alpha, 2, beta)*c0;
168  Real dtau_vu = elasticity_tensor(1, alpha, 0, beta)*c0;
169  Real dtau_vv = elasticity_tensor(1, alpha, 1, beta)*c0;
170  Real dtau_vw = elasticity_tensor(1, alpha, 2, beta)*c0;
171  Real dtau_wu = elasticity_tensor(2, alpha, 0, beta)*c0;
172  Real dtau_wv = elasticity_tensor(2, alpha, 1, beta)*c0;
173  Real dtau_ww = elasticity_tensor(2, alpha, 2, beta)*c0;
174 
175  Kuu(i,j) += dtau_uu*grad_phi[i][qp](alpha)*JxW[qp];
176  Kuv(i,j) += dtau_uv*grad_phi[i][qp](alpha)*JxW[qp];
177  Kuw(i,j) += dtau_uw*grad_phi[i][qp](alpha)*JxW[qp];
178  Kvu(i,j) += dtau_vu*grad_phi[i][qp](alpha)*JxW[qp];
179  Kvv(i,j) += dtau_vv*grad_phi[i][qp](alpha)*JxW[qp];
180  Kvw(i,j) += dtau_vw*grad_phi[i][qp](alpha)*JxW[qp];
181  Kwu(i,j) += dtau_wu*grad_phi[i][qp](alpha)*JxW[qp];
182  Kwv(i,j) += dtau_wv*grad_phi[i][qp](alpha)*JxW[qp];
183  Kww(i,j) += dtau_ww*grad_phi[i][qp](alpha)*JxW[qp];
184  }
185  }
186  }
187  }
188  }
189 
190  } // qp loop
191 
192  // If the Jacobian was requested, we computed it. Otherwise, we didn't.
193  return request_jacobian;
194 }
195 
196 bool ElasticitySystem::side_time_derivative (bool request_jacobian,
197  DiffContext & context)
198 {
199  FEMContext & c = cast_ref<FEMContext &>(context);
200 
201  // If we're on the correct side, apply the traction
202  if (c.has_side_boundary_id(BOUNDARY_ID_MAX_X))
203  {
204  // If we have an unsteady solver, then we need to extract the corresponding
205  // velocity variable. This allows us to use either a FirstOrderUnsteadySolver
206  // or a SecondOrderUnsteadySolver. That is, we get back the velocity variable
207  // index for FirstOrderUnsteadySolvers or, if it's a SecondOrderUnsteadySolver,
208  // this is actually just giving us back the same variable index.
209 
210  // If we only wanted to use a SecondOrderUnsteadySolver, then this
211  // step would be unnecessary and we would just
212  // populate the _u_var, etc. blocks of the residual and Jacobian.
213  unsigned int u_dot_var = this->get_second_order_dot_var(_u_var);
214  unsigned int v_dot_var = this->get_second_order_dot_var(_v_var);
215  unsigned int w_dot_var = this->get_second_order_dot_var(_w_var);
216 
217  FEBase * u_side_fe;
218  c.get_side_fe(_u_var, u_side_fe);
219 
220  // The number of local degrees of freedom in each variable
221  const unsigned int n_u_dofs = c.get_dof_indices(_u_var).size();
222 
223  DenseSubVector<Number> & Fu = c.get_elem_residual(u_dot_var);
224  DenseSubVector<Number> & Fv = c.get_elem_residual(v_dot_var);
225  DenseSubVector<Number> & Fw = c.get_elem_residual(w_dot_var);
226 
227  // Element Jacobian * quadrature weights for interior integration
228  const std::vector<Real> & JxW = u_side_fe->get_JxW();
229 
230  const std::vector<std::vector<Real>> & phi = u_side_fe->get_phi();
231 
232  unsigned int n_qpoints = c.get_side_qrule().n_points();
233 
234  Gradient traction(0.0, 0.0, -1.0);
235 
236  for (unsigned int qp=0; qp != n_qpoints; qp++)
237  {
238  for (unsigned int i=0; i != n_u_dofs; i++)
239  {
240  Fu(i) -= traction(0)*phi[i][qp]*JxW[qp];
241  Fv(i) -= traction(1)*phi[i][qp]*JxW[qp];
242  Fw(i) -= traction(2)*phi[i][qp]*JxW[qp];
243  }
244  }
245  }
246 
247  // If the Jacobian was requested, we computed it (in this case it's zero).
248  // Otherwise, we didn't.
249  return request_jacobian;
250 }
251 
252 bool ElasticitySystem::mass_residual(bool request_jacobian,
253  DiffContext & context)
254 {
255  FEMContext & c = cast_ref<FEMContext &>(context);
256 
257  // We need to extract the corresponding velocity variable.
258  // This allows us to use either a FirstOrderUnsteadySolver
259  // or a SecondOrderUnsteadySolver. That is, we get back the velocity variable
260  // index for FirstOrderUnsteadySolvers or, if it's a SecondOrderUnsteadySolver,
261  // this is actually just giving us back the same variable index.
262 
263  // If we only wanted to use a SecondOrderUnsteadySolver, then this
264  // step would be unnecessary and we would just
265  // populate the _u_var, etc. blocks of the residual and Jacobian.
266  unsigned int u_dot_var = this->get_second_order_dot_var(_u_var);
267  unsigned int v_dot_var = this->get_second_order_dot_var(_v_var);
268  unsigned int w_dot_var = this->get_second_order_dot_var(_w_var);
269 
270  FEBase * u_elem_fe;
271  c.get_element_fe(u_dot_var, u_elem_fe);
272 
273  // The number of local degrees of freedom in each variable
274  const unsigned int n_u_dofs = c.get_dof_indices(u_dot_var).size();
275 
276  // Element Jacobian * quadrature weights for interior integration
277  const std::vector<Real> & JxW = u_elem_fe->get_JxW();
278 
279  const std::vector<std::vector<Real>> & phi = u_elem_fe->get_phi();
280 
281  // Residuals that we're populating
285 
286  libMesh::DenseSubMatrix<libMesh::Number> & Kuu = c.get_elem_jacobian(u_dot_var, u_dot_var);
287  libMesh::DenseSubMatrix<libMesh::Number> & Kvv = c.get_elem_jacobian(v_dot_var, v_dot_var);
288  libMesh::DenseSubMatrix<libMesh::Number> & Kww = c.get_elem_jacobian(w_dot_var, w_dot_var);
289 
290  unsigned int n_qpoints = c.get_element_qrule().n_points();
291 
292  for (unsigned int qp=0; qp != n_qpoints; qp++)
293  {
294  // If we only cared about using FirstOrderUnsteadySolvers for time-stepping,
295  // then we could actually just use interior rate, but using interior_accel
296  // allows this assembly to function for both FirstOrderUnsteadySolvers
297  // and SecondOrderUnsteadySolvers
298  libMesh::Number u_ddot, v_ddot, w_ddot;
299  c.interior_accel(u_dot_var, qp, u_ddot);
300  c.interior_accel(v_dot_var, qp, v_ddot);
301  c.interior_accel(w_dot_var, qp, w_ddot);
302 
303  for (unsigned int i=0; i != n_u_dofs; i++)
304  {
305  Fu(i) += _rho*u_ddot*phi[i][qp]*JxW[qp];
306  Fv(i) += _rho*v_ddot*phi[i][qp]*JxW[qp];
307  Fw(i) += _rho*w_ddot*phi[i][qp]*JxW[qp];
308 
309  if (request_jacobian)
310  {
311  for (unsigned int j=0; j != n_u_dofs; j++)
312  {
313  libMesh::Real jac_term = _rho*phi[i][qp]*phi[j][qp]*JxW[qp];
314  jac_term *= context.get_elem_solution_accel_derivative();
315 
316  Kuu(i,j) += jac_term;
317  Kvv(i,j) += jac_term;
318  Kww(i,j) += jac_term;
319  }
320  }
321  }
322  }
323 
324  // If the Jacobian was requested, we computed it. Otherwise, we didn't.
325  return request_jacobian;
326 }
327 
328 Real ElasticitySystem::elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
329 {
330  // Hard code material parameters for the sake of simplicity
331  const Real poisson_ratio = 0.3;
332  const Real young_modulus = 1.0e2;
333 
334  // Define the Lame constants
335  const Real lambda_1 = (young_modulus*poisson_ratio)/((1.+poisson_ratio)*(1.-2.*poisson_ratio));
336  const Real lambda_2 = young_modulus/(2.*(1.+poisson_ratio));
337 
338  return
339  lambda_1 * kronecker_delta(i, j) * kronecker_delta(k, l) +
340  lambda_2 * (kronecker_delta(i, k) * kronecker_delta(j, l) + kronecker_delta(i, l) * kronecker_delta(j, k));
341 }
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
ConstFunction that simply returns 0.
Definition: zero_function.h:35
virtual void init_context(DiffContext &context)
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:772
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
Definition: assembly.C:38
virtual bool mass_residual(bool request_jacobian, DiffContext &context)
Subtracts a mass vector contribution on elem from elem_residual.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
The libMesh namespace provides an interface to certain functionality in the library.
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
Real get_elem_solution_accel_derivative() const
The derivative of the current elem_solution_accel w.r.t.
Definition: diff_context.h:437
Defines a dense subvector for use in finite element computations.
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
Defines a dense submatrix for use in Finite Element-type computations.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
void interior_accel(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1295
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
Real kronecker_delta(unsigned int i, unsigned int j)
Definition: assembly.C:32
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:381
virtual bool side_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on side of elem to elem_residual.
Real elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l)
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
bool has_side_boundary_id(boundary_id_type id) const
Reports if the boundary id is found on the current side.
Definition: fem_context.C:209
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Real get_elem_solution_derivative() const
The derivative of the current elem_solution w.r.t.
Definition: diff_context.h:419
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:845
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:299
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
unsigned int n_points() const
Definition: quadrature.h:113
sys get_dof_map()
This class forms the foundation from which generic finite elements may be derived.
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.