CArl
Code Arlequin / C++ implementation
assemble_functions_elasticity_3D.cpp
Go to the documentation of this file.
1 #include "assemble_functions_elasticity_3D.h"
2 
3 void set_homogeneous_physical_properties(libMesh::EquationSystems& es, std::string& physicalParamsFile)
4 {
5  const libMesh::Parallel::Communicator& SysComm = es.comm();
6  int rank = SysComm.rank();
7  int nodes = SysComm.size();
8 
9  // Read the random data info
10  double inputE;
11  double inputMu;
12 
13  if(rank == 0)
14  {
15  std::ifstream physicalParamsIFS(physicalParamsFile);
16  physicalParamsIFS >> inputE >> inputMu;
17  }
18 
19  if(nodes > 1)
20  {
21  SysComm.broadcast(inputE);
22  SysComm.broadcast(inputMu);
23  }
24 
25  // Mesh pointer
26  const libMesh::MeshBase& mesh = es.get_mesh();
27 
28  // Physical system and its "variables"
29  libMesh::ExplicitSystem& physical_param_system = es.get_system<libMesh::ExplicitSystem>("PhysicalConstants");
30  const libMesh::DofMap& physical_dof_map = physical_param_system.get_dof_map();
31 
32  unsigned int physical_consts[2];
33  physical_consts[0] = physical_param_system.variable_number ("E");
34  physical_consts[1] = physical_param_system.variable_number ("mu");
35 
36  std::vector<libMesh::dof_id_type> physical_dof_indices_var;
37  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
38  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
39 
40  for ( ; el != end_el; ++el)
41  {
42  const libMesh::Elem* elem = *el;
43 
44  // Young modulus, E
45  physical_dof_map.dof_indices(elem, physical_dof_indices_var, physical_consts[0]);
46  libMesh::dof_id_type dof_index = physical_dof_indices_var[0];
47 
48  if( (physical_param_system.solution->first_local_index() <= dof_index) &&
49  (dof_index < physical_param_system.solution->last_local_index()) )
50  {
51  physical_param_system.solution->set(dof_index, inputE);
52  }
53 
54  // Mu
55  physical_dof_map.dof_indices (elem, physical_dof_indices_var, physical_consts[1]);
56 
57  dof_index = physical_dof_indices_var[0];
58 
59  if( (physical_param_system.solution->first_local_index() <= dof_index) &&
60  (dof_index < physical_param_system.solution->last_local_index()) )
61  {
62  physical_param_system.solution->set(dof_index, inputMu);
63  }
64 
65  }
66 
67  physical_param_system.solution->close();
68  physical_param_system.update();
69 }
70 
71 void set_heterogeneous_physical_properties(libMesh::EquationSystems& es, std::string& physicalParamsFile)
72 {
73  const libMesh::Parallel::Communicator& SysComm = es.comm();
74  int rank = SysComm.rank();
75  int nodes = SysComm.size();
76 
77  // Read the random data info
78  std::vector<double> inputE;
79  std::vector<double> inputMu;
80  std::vector<int> inputIdx;
81 
82  int NbOfSubdomains = -1;
83 
84  double meanE = 0;
85  double meanMu = 0;
86 
87  if(rank == 0)
88  {
89  std::ifstream physicalParamsIFS(physicalParamsFile);
90  physicalParamsIFS >> NbOfSubdomains;
91  inputE.resize(NbOfSubdomains);
92  inputMu.resize(NbOfSubdomains);
93  inputIdx.resize(NbOfSubdomains);
94 
95  for(int iii = 0; iii < NbOfSubdomains; ++iii)
96  {
97  physicalParamsIFS >> inputE[iii];
98  physicalParamsIFS >> inputMu[iii];
99  physicalParamsIFS >> inputIdx[iii];
100 
101  meanE += inputE[iii];
102  meanMu += inputMu[iii];
103  }
104  meanE /= NbOfSubdomains;
105  meanMu /= NbOfSubdomains;
106  physicalParamsIFS.close();
107  }
108 
109  if(nodes > 1)
110  {
111  SysComm.broadcast(NbOfSubdomains);
112  SysComm.broadcast(meanE);
113  SysComm.broadcast(meanMu);
114 
115  if(rank != 0)
116  {
117  inputE.resize(NbOfSubdomains);
118  inputMu.resize(NbOfSubdomains);
119  inputIdx.resize(NbOfSubdomains);
120  }
121  SysComm.broadcast(inputE);
122  SysComm.broadcast(inputMu);
123  SysComm.broadcast(inputIdx);
124  }
125 
126  // Mesh pointer
127  const libMesh::MeshBase& mesh = es.get_mesh();
128 
129  // Physical system and its "variables"
130  libMesh::ExplicitSystem& physical_param_system = es.get_system<libMesh::ExplicitSystem>("PhysicalConstants");
131  const libMesh::DofMap& physical_dof_map = physical_param_system.get_dof_map();
132 
133  unsigned int physical_consts[2];
134  physical_consts[0] = physical_param_system.variable_number ("E");
135  physical_consts[1] = physical_param_system.variable_number ("mu");
136 
137  std::vector<libMesh::dof_id_type> physical_dof_indices_var;
138  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
139  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
140 
141  int currentSubdomain = -1;
142 
143  for ( ; el != end_el; ++el)
144  {
145  const libMesh::Elem* elem = *el;
146 
147  currentSubdomain = elem->subdomain_id()-1;
148 
149  // Young modulus, E
150  physical_dof_map.dof_indices(elem, physical_dof_indices_var, physical_consts[0]);
151  libMesh::dof_id_type dof_index = physical_dof_indices_var[0];
152 
153  if( (physical_param_system.solution->first_local_index() <= dof_index) &&
154  (dof_index < physical_param_system.solution->last_local_index()) )
155  {
156  physical_param_system.solution->set(dof_index, inputE[currentSubdomain]);
157  }
158 
159  // Mu
160  physical_dof_map.dof_indices (elem, physical_dof_indices_var, physical_consts[1]);
161 
162  dof_index = physical_dof_indices_var[0];
163 
164  if( (physical_param_system.solution->first_local_index() <= dof_index) &&
165  (dof_index < physical_param_system.solution->last_local_index()) )
166  {
167  physical_param_system.solution->set(dof_index, inputMu[currentSubdomain]);
168  }
169  }
170 
171  physical_param_system.solution->close();
172  physical_param_system.update();
173 }
174 
175 libMesh::Real eval_lambda_1(libMesh::Real E, libMesh::Real mu)
176 {
177  return mu*(E - 2*mu)/(3*mu-E);
178 }
179 
180 libMesh::Real eval_elasticity_tensor(unsigned int i,
181  unsigned int j,
182  unsigned int k,
183  unsigned int l,
184  libMesh::Number E,
185  libMesh::Number mu)
186 {
187  const libMesh::Real lambda_1 = eval_lambda_1(E,mu);
188  const libMesh::Real lambda_2 = mu;
189 
190  return lambda_1 * kronecker_delta(i,j) * kronecker_delta(k,l)
191  + lambda_2 * (kronecker_delta(i,k) * kronecker_delta(j,l)
192  + kronecker_delta(i,l) * kronecker_delta(j,k));
193 }
194 
195 void Update_SubK_isotropic( libMesh::DenseSubMatrix<libMesh::Number>& SubK,
196  unsigned int qp,
197  unsigned int C_i,
198  unsigned int C_k,
199  const std::vector<std::vector<libMesh::RealGradient> >& dphi,
200  const unsigned int n_components,
201  const unsigned int n_u_dofs,
202  const std::vector<libMesh::Real>& JxW,
203  libMesh::Number E,
204  libMesh::Number mu,
205  double cte
206  )
207 {
208  for (unsigned int iii=0; iii<n_u_dofs; iii++)
209  {
210  for (unsigned int jjj=0; jjj<n_u_dofs; jjj++)
211  {
212  for(unsigned int C_j=0; C_j<n_components; C_j++)
213  {
214  for(unsigned int C_l=0; C_l<n_components; C_l++)
215  {
216  SubK(iii,jjj) += cte * JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l,E,mu) * dphi[iii][qp](C_j)*dphi[jjj][qp](C_l));
217  }
218  }
219  }
220  }
221 };
222 
223 void assemble_elasticity_with_weight( libMesh::EquationSystems& es,
224  const std::string& system_name,
225  weight_parameter_function& weight_mask,
226  WeightFunctionSystemType system_type)
227 {
228  libmesh_assert_equal_to(system_name, "Elasticity");
229 
230  libMesh::PerfLog perf_log ("Matrix Assembly (Homogeneous elasticity)",MASTER_bPerfLog_assemble_fem);
231 
232  perf_log.push("Preamble");
233 
234  const libMesh::MeshBase& mesh = es.get_mesh();
235 
236  const unsigned int dim = mesh.mesh_dimension();
237 
238  // - Set up physical properties system ------------------------------------
239  libMesh::ExplicitSystem& physical_param_system = es.get_system<libMesh::ExplicitSystem>("PhysicalConstants");
240  const unsigned int young_var = physical_param_system.variable_number ("E");
241  const unsigned int mu_var = physical_param_system.variable_number ("mu");
242 
243  const libMesh::DofMap& physical_dof_map = physical_param_system.get_dof_map();
244  std::vector<libMesh::dof_id_type> physical_dof_indices_var;
245 
246  // The DoF and values of the physical system
247  const libMesh::Elem* phys_elem = *(mesh.active_local_elements_begin());
248  physical_dof_map.dof_indices(phys_elem, physical_dof_indices_var, young_var);
249  libMesh::Number localE = physical_param_system.current_solution(physical_dof_indices_var[0]);
250 
251  physical_dof_map.dof_indices(phys_elem, physical_dof_indices_var, mu_var);
252  libMesh::Number localMu = physical_param_system.current_solution(physical_dof_indices_var[0]);
253 
254  // - Set up elasticity system ---------------------------------------------
255  libMesh::LinearImplicitSystem& system = es.get_system<libMesh::LinearImplicitSystem>("Elasticity");
256 
257  const unsigned int n_components = 3;
258  const unsigned int u_var = system.variable_number ("u");
259  const unsigned int v_var = system.variable_number ("v");
260  const unsigned int w_var = system.variable_number ("w");
261 
262  const libMesh::DofMap& dof_map = system.get_dof_map();
263  libMesh::FEType fe_type = dof_map.variable_type(u_var);
264 
265  // Set up pointers to FEBase's of dimension dim and FE type fe_type
266  // -> 3D elements
267  libMesh::UniquePtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
268  libMesh::QGauss qrule (dim, fe_type.default_quadrature_order());
269  fe->attach_quadrature_rule (&qrule);
270 
271  // -> Faces
272  libMesh::UniquePtr<libMesh::FEBase> fe_face (libMesh::FEBase::build(dim, fe_type));
273  libMesh::QGauss qface(dim-1, fe_type.default_quadrature_order());
274  fe_face->attach_quadrature_rule (&qface);
275 
276  // Jacobian
277  const std::vector<libMesh::Real>& JxW = fe->get_JxW();
278 
279  // Shape functions derivatives
280  const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
281 
282  // Quadrature points
283  const std::vector<libMesh::Point>& qp_points = fe->get_xyz();
284 
285  // Weights for the Arlequin method
286  double weight_alpha = 1;
287 
288  libMesh::DenseMatrix<libMesh::Number> Ke;
289  libMesh::DenseVector<libMesh::Number> Fe;
290 
291  libMesh::DenseSubMatrix<libMesh::Number>
292  Kuu(Ke), Kuv(Ke), Kuw(Ke),
293  Kvu(Ke), Kvv(Ke), Kvw(Ke),
294  Kwu(Ke), Kwv(Ke), Kww(Ke);
295 
296  libMesh::DenseSubVector<libMesh::Number>
297  Fu(Fe),
298  Fv(Fe),
299  Fw(Fe);
300 
301  std::vector<libMesh::dof_id_type> dof_indices;
302  std::vector<libMesh::dof_id_type> dof_indices_u;
303  std::vector<libMesh::dof_id_type> dof_indices_v;
304  std::vector<libMesh::dof_id_type> dof_indices_w;
305 
306  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
307  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
308 
309  perf_log.pop("Preamble");
310 
311  // For each element
312  for ( ; el != end_el; ++el)
313  {
314  // Get its pointer
315  const libMesh::Elem* elem = *el;
316 
317  perf_log.push("Define DoF");
318  // The total DoF indices, and those associated to each variable
319  dof_map.dof_indices (elem, dof_indices);
320  dof_map.dof_indices (elem, dof_indices_u, u_var);
321  dof_map.dof_indices (elem, dof_indices_v, v_var);
322  dof_map.dof_indices (elem, dof_indices_w, w_var);
323 
324  const unsigned int n_dofs = dof_indices.size();
325  const unsigned int n_u_dofs = dof_indices_u.size();
326  const unsigned int n_v_dofs = dof_indices_v.size();
327  const unsigned int n_w_dofs = dof_indices_w.size();
328 
329  perf_log.pop("Define DoF");
330 
331  // Restart the FE to the "geometry" of the element
332  // -> Determines quadrature points, shape functions ...
333  // !!! User can change the points to be used (can use other mesh's points
334  // instead of the quadrature points)
335  fe->reinit (elem);
336 
337  perf_log.push("Matrix manipulations");
338  Ke.resize (n_dofs, n_dofs);
339  Fe.resize (n_dofs);
340 
341  // Set the positions of the sub-matrices
342  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
343  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
344  Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
345 
346  Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
347  Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
348  Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
349 
350  Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
351  Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
352  Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
353 
354  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
355  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
356  Fw.reposition (w_var*n_u_dofs, n_w_dofs);
357  perf_log.push("Matrix manipulations");
358 
359  // For each quadrature point determinate the sub-matrices elements
360  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
361  {
362 
363  perf_log.push("Rigidity","Matrix element calculations");
364  weight_alpha = weight_mask.get_alpha(qp_points[qp],system_type);
365 
366  // Internal tension
367  Update_SubK_isotropic(Kuu, qp, 0, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
368  Update_SubK_isotropic(Kuv, qp, 0, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
369  Update_SubK_isotropic(Kuw, qp, 0, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
370 
371  Update_SubK_isotropic(Kvu, qp, 1, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
372  Update_SubK_isotropic(Kvv, qp, 1, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
373  Update_SubK_isotropic(Kvw, qp, 1, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
374 
375  Update_SubK_isotropic(Kwu, qp, 2, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
376  Update_SubK_isotropic(Kwv, qp, 2, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
377  Update_SubK_isotropic(Kww, qp, 2, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
378 
379  perf_log.pop("Rigidity","Matrix element calculations");
380  }
381 
382  // Apply constraints
383  perf_log.push("Constraints","Matrix element calculations");
384  dof_map.heterogenously_constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
385  perf_log.pop("Constraints","Matrix element calculations");
386 
387  perf_log.push("Adding elements");
388  system.matrix->add_matrix (Ke, dof_indices);
389  system.rhs->add_vector (Fe, dof_indices);
390  perf_log.pop("Adding elements");
391  }
392 
393  system.matrix->close();
394  system.rhs->close();
395 }
396 
397 void assemble_elasticity_with_weight_and_traction(libMesh::EquationSystems& es,
398  const std::string& system_name,
399  weight_parameter_function& weight_mask,
400  WeightFunctionSystemType system_type,
401  int traction_boundary_id,
402  std::vector<double> traction_density)
403 
404 {
405  libmesh_assert_equal_to(system_name, "Elasticity");
406 
407  libMesh::PerfLog perf_log ("Matrix Assembly (Homogeneous elasticity)",MASTER_bPerfLog_assemble_fem);
408 
409  perf_log.push("Preamble");
410 
411  const libMesh::MeshBase& mesh = es.get_mesh();
412 
413  const unsigned int dim = mesh.mesh_dimension();
414 
415  // - Set up physical properties system ------------------------------------
416  libMesh::ExplicitSystem& physical_param_system = es.get_system<libMesh::ExplicitSystem>("PhysicalConstants");
417  const unsigned int young_var = physical_param_system.variable_number ("E");
418  const unsigned int mu_var = physical_param_system.variable_number ("mu");
419 
420  const libMesh::DofMap& physical_dof_map = physical_param_system.get_dof_map();
421  std::vector<libMesh::dof_id_type> physical_dof_indices_var;
422 
423  // The DoF and values of the physical system
424  const libMesh::Elem* phys_elem = *(mesh.active_local_elements_begin());
425  physical_dof_map.dof_indices(phys_elem, physical_dof_indices_var, young_var);
426  libMesh::Number localE = physical_param_system.current_solution(physical_dof_indices_var[0]);
427 
428  physical_dof_map.dof_indices(phys_elem, physical_dof_indices_var, mu_var);
429  libMesh::Number localMu = physical_param_system.current_solution(physical_dof_indices_var[0]);
430 
431  // - Set up elasticity system ---------------------------------------------
432  libMesh::LinearImplicitSystem& system = es.get_system<libMesh::LinearImplicitSystem>("Elasticity");
433 
434  const unsigned int n_components = 3;
435  const unsigned int u_var = system.variable_number ("u");
436  const unsigned int v_var = system.variable_number ("v");
437  const unsigned int w_var = system.variable_number ("w");
438 
439  const libMesh::DofMap& dof_map = system.get_dof_map();
440  libMesh::FEType fe_type = dof_map.variable_type(u_var);
441 
442  // Set up pointers to FEBase's of dimension dim and FE type fe_type
443  // -> 3D elements
444  libMesh::UniquePtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
445  libMesh::QGauss qrule (dim, fe_type.default_quadrature_order());
446  fe->attach_quadrature_rule (&qrule);
447 
448  // -> Faces
449  libMesh::UniquePtr<libMesh::FEBase> fe_face (libMesh::FEBase::build(dim, fe_type));
450  libMesh::QGauss qface(dim-1, fe_type.default_quadrature_order());
451  fe_face->attach_quadrature_rule (&qface);
452 
453  // -> Face traction
454  libMesh::DenseVector<libMesh::Number> g_vec(3);
455  g_vec(0) = traction_density[0];
456  g_vec(1) = traction_density[1];
457  g_vec(2) = traction_density[2];
458 
459  // Jacobian
460  const std::vector<libMesh::Real>& JxW = fe->get_JxW();
461 
462  // Shape functions derivatives
463  const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
464 
465  // Quadrature points
466  const std::vector<libMesh::Point>& qp_points = fe->get_xyz();
467 
468  // Face Jacobian
469  const std::vector<libMesh::Real> & JxW_face = fe_face->get_JxW();
470 
471  // Face shape functions
472  const std::vector<std::vector<libMesh::Real> > & phi_face = fe_face->get_phi();
473 
474  // Face quadrature points
475  const std::vector<libMesh::Point>& qp_face_points = fe_face->get_xyz();
476 
477 
478  // Weights for the Arlequin method
479  double weight_alpha = 1;
480 
481  libMesh::DenseMatrix<libMesh::Number> Ke;
482  libMesh::DenseVector<libMesh::Number> Fe;
483 
484  libMesh::DenseSubMatrix<libMesh::Number>
485  Kuu(Ke), Kuv(Ke), Kuw(Ke),
486  Kvu(Ke), Kvv(Ke), Kvw(Ke),
487  Kwu(Ke), Kwv(Ke), Kww(Ke);
488 
489  libMesh::DenseSubVector<libMesh::Number>
490  Fu(Fe),
491  Fv(Fe),
492  Fw(Fe);
493 
494  std::vector<libMesh::dof_id_type> dof_indices;
495  std::vector<libMesh::dof_id_type> dof_indices_u;
496  std::vector<libMesh::dof_id_type> dof_indices_v;
497  std::vector<libMesh::dof_id_type> dof_indices_w;
498 
499  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
500  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
501 
502  perf_log.pop("Preamble");
503 
504  // For each element
505  for ( ; el != end_el; ++el)
506  {
507  // Get its pointer
508  const libMesh::Elem* elem = *el;
509 
510  perf_log.push("Define DoF");
511  // The total DoF indices, and those associated to each variable
512  dof_map.dof_indices (elem, dof_indices);
513  dof_map.dof_indices (elem, dof_indices_u, u_var);
514  dof_map.dof_indices (elem, dof_indices_v, v_var);
515  dof_map.dof_indices (elem, dof_indices_w, w_var);
516 
517  const unsigned int n_dofs = dof_indices.size();
518  const unsigned int n_u_dofs = dof_indices_u.size();
519  const unsigned int n_v_dofs = dof_indices_v.size();
520  const unsigned int n_w_dofs = dof_indices_w.size();
521 
522  perf_log.pop("Define DoF");
523 
524  // Restart the FE to the "geometry" of the element
525  // -> Determines quadrature points, shape functions ...
526  // !!! User can change the points to be used (can use other mesh's points
527  // instead of the quadrature points)
528  fe->reinit (elem);
529 
530  perf_log.push("Matrix manipulations");
531  Ke.resize (n_dofs, n_dofs);
532  Fe.resize (n_dofs);
533 
534  // Set the positions of the sub-matrices
535  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
536  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
537  Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
538 
539  Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
540  Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
541  Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
542 
543  Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
544  Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
545  Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
546 
547  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
548  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
549  Fw.reposition (w_var*n_u_dofs, n_w_dofs);
550  perf_log.push("Matrix manipulations");
551 
552  // For each quadrature point determinate the sub-matrices elements
553  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
554  {
555 
556  perf_log.push("Rigidity","Matrix element calculations");
557  weight_alpha = weight_mask.get_alpha(qp_points[qp],system_type);
558 
559  // Internal tension
560  Update_SubK_isotropic(Kuu, qp, 0, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
561  Update_SubK_isotropic(Kuv, qp, 0, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
562  Update_SubK_isotropic(Kuw, qp, 0, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
563 
564  Update_SubK_isotropic(Kvu, qp, 1, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
565  Update_SubK_isotropic(Kvv, qp, 1, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
566  Update_SubK_isotropic(Kvw, qp, 1, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
567 
568  Update_SubK_isotropic(Kwu, qp, 2, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
569  Update_SubK_isotropic(Kwv, qp, 2, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
570  Update_SubK_isotropic(Kww, qp, 2, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
571 
572  perf_log.pop("Rigidity","Matrix element calculations");
573  }
574 
575  // Assemble the traction
576  for (unsigned int side=0; side<elem->n_sides(); side++)
577  {
578  if (elem->neighbor_ptr(side) == libmesh_nullptr)
579  {
580  fe_face->reinit(elem, side);
581 
582  // Apply a traction
583  for (unsigned int qp=0; qp<qface.n_points(); qp++)
584  {
585  weight_alpha = weight_mask.get_alpha(qp_face_points[qp],system_type);
586 
587  if (mesh.get_boundary_info().has_boundary_id(elem, side, traction_boundary_id))
588  {
589  for (unsigned int dof_iii=0; dof_iii<n_u_dofs; dof_iii++)
590  {
591  Fu(dof_iii) += weight_alpha * JxW_face[qp] * (g_vec(0) * phi_face[dof_iii][qp]);
592  Fv(dof_iii) += weight_alpha * JxW_face[qp] * (g_vec(1) * phi_face[dof_iii][qp]);
593  Fw(dof_iii) += weight_alpha * JxW_face[qp] * (g_vec(2) * phi_face[dof_iii][qp]);
594  }
595  }
596  }
597  }
598  }
599 
600  // Apply constraints
601  perf_log.push("Constraints","Matrix element calculations");
602  dof_map.heterogenously_constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
603  perf_log.pop("Constraints","Matrix element calculations");
604 
605  perf_log.push("Adding elements");
606  system.matrix->add_matrix (Ke, dof_indices);
607  system.rhs->add_vector (Fe, dof_indices);
608  perf_log.pop("Adding elements");
609  }
610 
611  system.matrix->close();
612  system.rhs->close();
613 };
614 
615 void assemble_elasticity_heterogeneous_with_weight( libMesh::EquationSystems& es,
616  const std::string& system_name,
617  weight_parameter_function& weight_mask,
618  WeightFunctionSystemType system_type)
619 {
620  libmesh_assert_equal_to (system_name, "Elasticity");
621 
622  libMesh::PerfLog perf_log ("Matrix Assembly (Heterogeneous elasticity)",MASTER_bPerfLog_assemble_fem);
623 
624  perf_log.push("Preamble");
625  // Set up mesh
626  const libMesh::MeshBase& mesh = es.get_mesh();
627 
628  const unsigned int dim = mesh.mesh_dimension();
629 
630  // - Set up physical properties system ------------------------------------
631  libMesh::Number localE = -1;
632  libMesh::Number localMu = -1;
633 
634  libMesh::ExplicitSystem& physical_param_system = es.get_system<libMesh::ExplicitSystem>("PhysicalConstants");
635  const unsigned int young_var = physical_param_system.variable_number ("E");
636  const unsigned int mu_var = physical_param_system.variable_number ("mu");
637 
638  const libMesh::DofMap& physical_dof_map = physical_param_system.get_dof_map();
639  std::vector<libMesh::dof_id_type> physical_dof_indices_var;
640 
641  // - Set up elasticity system ---------------------------------------------
642  libMesh::LinearImplicitSystem& system = es.get_system<libMesh::LinearImplicitSystem>("Elasticity");
643 
644  const unsigned int n_components = 3;
645  const unsigned int u_var = system.variable_number ("u");
646  const unsigned int v_var = system.variable_number ("v");
647  const unsigned int w_var = system.variable_number ("w");
648 
649  const libMesh::DofMap& dof_map = system.get_dof_map();
650  libMesh::FEType fe_type = dof_map.variable_type(u_var);
651 
652  // Set up pointers to FEBase's of dimension dim and FE type fe_type
653  // -> 3D elements
654  libMesh::UniquePtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
655  libMesh::QGauss qrule (dim, fe_type.default_quadrature_order());
656  fe->attach_quadrature_rule (&qrule);
657 
658  // -> Faces
659  libMesh::UniquePtr<libMesh::FEBase> fe_face (libMesh::FEBase::build(dim, fe_type));
660  libMesh::QGauss qface(dim-1, fe_type.default_quadrature_order());
661  fe_face->attach_quadrature_rule (&qface);
662 
663  // Jacobian
664  const std::vector<libMesh::Real>& JxW = fe->get_JxW();
665 
666  // Shape functions derivatives
667  const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
668 
669  // Quadrature points
670  const std::vector<libMesh::Point>& qp_points = fe->get_xyz();
671 
672  // Weights for the Arlequin method
673  double weight_alpha = 1;
674 
675  libMesh::DenseMatrix<libMesh::Number> Ke;
676  libMesh::DenseVector<libMesh::Number> Fe;
677 
678  libMesh::DenseSubMatrix<libMesh::Number>
679  Kuu(Ke), Kuv(Ke), Kuw(Ke),
680  Kvu(Ke), Kvv(Ke), Kvw(Ke),
681  Kwu(Ke), Kwv(Ke), Kww(Ke);
682 
683  libMesh::DenseSubVector<libMesh::Number>
684  Fu(Fe),
685  Fv(Fe),
686  Fw(Fe);
687 
688  std::vector<libMesh::dof_id_type> dof_indices;
689  std::vector<libMesh::dof_id_type> dof_indices_u;
690  std::vector<libMesh::dof_id_type> dof_indices_v;
691  std::vector<libMesh::dof_id_type> dof_indices_w;
692 
693  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
694  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
695 
696  perf_log.pop("Preamble");
697 
698  // For each element
699  for ( ; el != end_el; ++el)
700  {
701  // Get its pointer
702  const libMesh::Elem* elem = *el;
703 
704  perf_log.push("Define DoF");
705  // The total DoF indices, and those associated to each variable
706  dof_map.dof_indices (elem, dof_indices);
707  dof_map.dof_indices (elem, dof_indices_u, u_var);
708  dof_map.dof_indices (elem, dof_indices_v, v_var);
709  dof_map.dof_indices (elem, dof_indices_w, w_var);
710 
711  const unsigned int n_dofs = dof_indices.size();
712  const unsigned int n_u_dofs = dof_indices_u.size();
713  const unsigned int n_v_dofs = dof_indices_v.size();
714  const unsigned int n_w_dofs = dof_indices_w.size();
715 
716  perf_log.pop("Define DoF");
717 
718  perf_log.push("Define physical params");
719  // The DoF and values of the physical system
720  physical_dof_map.dof_indices(elem, physical_dof_indices_var, young_var);
721  localE = physical_param_system.current_solution(physical_dof_indices_var[0]);
722 
723  physical_dof_map.dof_indices(elem, physical_dof_indices_var, mu_var);
724  localMu = physical_param_system.current_solution(physical_dof_indices_var[0]);
725  perf_log.pop("Define physical params");
726 
727  // Restart the FE to the "geometry" of the element
728  // -> Determines quadrature points, shape functions ...
729  // !!! User can change the points to be used (can use other mesh's points
730  // instead of the quadrature points)
731  fe->reinit (elem);
732 
733  perf_log.push("Matrix manipulations");
734  Ke.resize (n_dofs, n_dofs);
735  Fe.resize (n_dofs);
736 
737  // Set the positions of the sub-matrices
738  Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
739  Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
740  Kuw.reposition (u_var*n_u_dofs, w_var*n_u_dofs, n_u_dofs, n_w_dofs);
741 
742  Kvu.reposition (v_var*n_u_dofs, u_var*n_u_dofs, n_v_dofs, n_u_dofs);
743  Kvv.reposition (v_var*n_u_dofs, v_var*n_u_dofs, n_v_dofs, n_v_dofs);
744  Kvw.reposition (v_var*n_u_dofs, w_var*n_u_dofs, n_v_dofs, n_w_dofs);
745 
746  Kwu.reposition (w_var*n_u_dofs, u_var*n_u_dofs, n_w_dofs, n_u_dofs);
747  Kwv.reposition (w_var*n_u_dofs, v_var*n_u_dofs, n_w_dofs, n_v_dofs);
748  Kww.reposition (w_var*n_u_dofs, w_var*n_u_dofs, n_w_dofs, n_w_dofs);
749 
750  Fu.reposition (u_var*n_u_dofs, n_u_dofs);
751  Fv.reposition (v_var*n_u_dofs, n_v_dofs);
752  Fw.reposition (w_var*n_u_dofs, n_w_dofs);
753  perf_log.push("Matrix manipulations");
754 
755  perf_log.push("Matrix element calculations");
756  // For each quadrature point determinate the sub-matrices elements
757  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
758  {
759 
760  perf_log.push("Rigidity","Matrix element calculations");
761  weight_alpha = weight_mask.get_alpha(qp_points[qp],system_type);
762 
763  // Internal tension
764  Update_SubK_isotropic(Kuu, qp, 0, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
765  Update_SubK_isotropic(Kuv, qp, 0, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
766  Update_SubK_isotropic(Kuw, qp, 0, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
767 
768  Update_SubK_isotropic(Kvu, qp, 1, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
769  Update_SubK_isotropic(Kvv, qp, 1, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
770  Update_SubK_isotropic(Kvw, qp, 1, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
771 
772  Update_SubK_isotropic(Kwu, qp, 2, 0, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
773  Update_SubK_isotropic(Kwv, qp, 2, 1, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
774  Update_SubK_isotropic(Kww, qp, 2, 2, dphi, n_components, n_u_dofs, JxW, localE, localMu, weight_alpha);
775 
776  perf_log.pop("Rigidity","Matrix element calculations");
777  }
778 
779  // Apply constraints
780  perf_log.push("Constraints","Matrix element calculations");
781  dof_map.heterogenously_constrain_element_matrix_and_vector(Ke, Fe, dof_indices);
782  perf_log.pop("Constraints","Matrix element calculations");
783 
784  perf_log.push("Adding elements");
785  system.matrix->add_matrix (Ke, dof_indices);
786  system.rhs->add_vector (Fe, dof_indices);
787  perf_log.pop("Adding elements");
788  }
789 
790  system.matrix->close();
791  system.rhs->close();
792 };
793 
794 void compute_stresses(libMesh::EquationSystems& es)
795 {
796  // Get mesh pointer
797  const libMesh::MeshBase& mesh = es.get_mesh();
798 
799  const unsigned int dim = mesh.mesh_dimension();
800 
801  // Get the Elasticity system
802  libMesh::LinearImplicitSystem& elasticity_system = es.get_system<libMesh::LinearImplicitSystem>("Elasticity");
803 
804  // Set the displacement variables
805  unsigned int displacement_vars[3];
806  displacement_vars[0] = elasticity_system.variable_number ("u");
807  displacement_vars[1] = elasticity_system.variable_number ("v");
808  displacement_vars[2] = elasticity_system.variable_number ("w");
809  const unsigned int u_var = elasticity_system.variable_number ("u");
810 
811  // - Set up physical properties system ------------------------------------
812  libMesh::ExplicitSystem& physical_param_system = es.get_system<libMesh::ExplicitSystem>("PhysicalConstants");
813  const unsigned int young_var = physical_param_system.variable_number ("E");
814  const unsigned int mu_var = physical_param_system.variable_number ("mu");
815 
816  const libMesh::DofMap& physical_dof_map = physical_param_system.get_dof_map();
817  std::vector<libMesh::dof_id_type> physical_dof_indices_var;
818 
819  // The DoF and values of the physical system
820  const libMesh::Elem* phys_elem = *(mesh.active_local_elements_begin());
821  physical_dof_map.dof_indices(phys_elem, physical_dof_indices_var, young_var);
822  libMesh::Number localE = physical_param_system.current_solution(physical_dof_indices_var[0]);
823 
824  physical_dof_map.dof_indices(phys_elem, physical_dof_indices_var, mu_var);
825  libMesh::Number localMu = physical_param_system.current_solution(physical_dof_indices_var[0]);
826 
827  // - Degrees of freedom ---------------------------------------------------
828  // Get the DoF of the Elasticity system
829  const libMesh::DofMap& elasticity_dof_map = elasticity_system.get_dof_map();
830 
831  // Get the type of FE used for "u_var" (same as others)
832  libMesh::FEType fe_type = elasticity_dof_map.variable_type(u_var);
833 
834  // Build a FE to be used here and attach a default quadrature
835  libMesh::UniquePtr<libMesh::FEBase> fe (libMesh::FEBase::build(dim, fe_type));
836  libMesh::QGauss qrule (dim, fe_type.default_quadrature_order());
837  fe->attach_quadrature_rule (&qrule);
838 
839  // - Stress variables
840 
841  // Jacobian vector and d_phi
842  const std::vector<libMesh::Real>& JxW = fe->get_JxW();
843  const std::vector<std::vector<libMesh::RealGradient> >& dphi = fe->get_dphi();
844 
845  // Get the stress system and its DoF map
846  libMesh::ExplicitSystem& stress_system = es.get_system<libMesh::ExplicitSystem>("StressSystem");
847  const libMesh::DofMap& stress_dof_map = stress_system.get_dof_map();
848 
849  // The 9 stress variables + von Mises
850  unsigned int sigma_vars[3][3];
851  sigma_vars[0][0] = stress_system.variable_number ("sigma_00");
852  sigma_vars[0][1] = stress_system.variable_number ("sigma_01");
853  sigma_vars[0][2] = stress_system.variable_number ("sigma_02");
854  sigma_vars[1][0] = stress_system.variable_number ("sigma_10");
855  sigma_vars[1][1] = stress_system.variable_number ("sigma_11");
856  sigma_vars[1][2] = stress_system.variable_number ("sigma_12");
857  sigma_vars[2][0] = stress_system.variable_number ("sigma_20");
858  sigma_vars[2][1] = stress_system.variable_number ("sigma_21");
859  sigma_vars[2][2] = stress_system.variable_number ("sigma_22");
860  unsigned int vonMises_var = stress_system.variable_number ("vonMises");
861 
862  // Define a table of DoF indices for each elasticity variable
863  std::vector< std::vector<libMesh::dof_id_type> > elasticity_dof_indices_var(elasticity_system.n_vars());
864 
865  // Define a vector for the stress DoF indices
866  std::vector<libMesh::dof_id_type> stress_dof_indices_var;
867 
868  // "Local" matrix
869  libMesh::DenseMatrix<libMesh::Number> elem_sigma;
870 
871  // Iterators over the local elements
872  libMesh::MeshBase::const_element_iterator el = mesh.active_local_elements_begin();
873  const libMesh::MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
874 
875  for ( ; el != end_el; ++el)
876  {
877  const libMesh::Elem* elem = *el;
878 
879  /* For each variable identified by "displacement_vars[var]", get
880  * the GLOBAL indices from the elasticity DoF map and save it in the
881  * vector "elasticity_dof_indices_var[var]"
882  */
883  for(unsigned int var=0; var<3; var++)
884  {
885  elasticity_dof_map.dof_indices (elem, elasticity_dof_indices_var[var], displacement_vars[var]);
886  }
887 
888  physical_dof_map.dof_indices(elem, physical_dof_indices_var, young_var);
889  localE = physical_param_system.current_solution(physical_dof_indices_var[0]);
890 
891  physical_dof_map.dof_indices(elem, physical_dof_indices_var, mu_var);
892  localMu = physical_param_system.current_solution(physical_dof_indices_var[0]);
893 
894  // Restart the element properties and the local stress contribution matrix
895  fe->reinit (elem);
896  elem_sigma.resize(3,3);
897 
898  // For each quadrature point
899  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
900  {
901  for(unsigned int C_i=0; C_i<3; C_i++)
902  {
903  for(unsigned int C_j=0; C_j<3; C_j++)
904  {
905  for(unsigned int C_k=0; C_k<3; C_k++)
906  {
907  // Number of DoFs of the elasticity
908  const unsigned int n_var_dofs = elasticity_dof_indices_var[C_k].size();
909 
910  libMesh::Gradient displacement_gradient_k;
911  for(unsigned int l=0; l<n_var_dofs; l++)
912  {
913  /*
914  * For each DoF of the elasticity problem, "l",
915  * get the gradient of the shape functions at the
916  * quadrature point "qp", and add it to the total
917  * displacement gradient, scaled by the elasticity's
918  * current solution for the displacement variable
919  * "C_k"
920  */
921  displacement_gradient_k.add_scaled(dphi[l][qp], elasticity_system.current_solution(elasticity_dof_indices_var[C_k][l]));
922  }
923 
924  for(unsigned int C_l=0; C_l<3; C_l++)
925  {
926  // S_{i,j} = C_{i,j,k,l} * E_{k,l}
927  elem_sigma(C_i,C_j) += JxW[qp]*(eval_elasticity_tensor(C_i,C_j,C_k,C_l,localE,localMu) * displacement_gradient_k(C_l) );
928  }
929  }
930  }
931  }
932  }
933 
934  // Rescale the matrix by the elem. volume
935  elem_sigma.scale(1./elem->volume());
936 
937  for(unsigned int i=0; i<3; i++)
938  {
939  for(unsigned int j=0; j<3; j++)
940  {
941  // For each variable of the stress system, get the DoF map
942  stress_dof_map.dof_indices(elem, stress_dof_indices_var, sigma_vars[i][j]);
943 
944  // Get the first index (CONSTANT MONOMIAL basis functions, hence
945  // only one element
946  libMesh::dof_id_type dof_index = stress_dof_indices_var[0];
947 
948  // To be sure, test if the DoF index is inside this processor
949  if( (stress_system.solution->first_local_index() <= dof_index) &&
950  (dof_index < stress_system.solution->last_local_index()) )
951  {
952  stress_system.solution->set(dof_index, elem_sigma(i,j));
953  }
954  }
955  }
956 
957  // Calculate von Mises
958  libMesh::Number vonMises_value = std::sqrt( 0.5*( pow(elem_sigma(0,0) - elem_sigma(1,1),2.) +
959  pow(elem_sigma(1,1) - elem_sigma(2,2),2.) +
960  pow(elem_sigma(2,2) - elem_sigma(0,0),2.) +
961  6.*(pow(elem_sigma(0,1),2.) + pow(elem_sigma(1,2),2.) + pow(elem_sigma(2,0),2.))
962  ) );
963 
964  // Get the DoF map
965  stress_dof_map.dof_indices (elem, stress_dof_indices_var, vonMises_var);
966 
967  // Get the first index (CONSTANT MONOMIAL basis functions, hence
968  // only one element
969  libMesh::dof_id_type dof_index = stress_dof_indices_var[0];
970 
971  // To be sure, test if the DoF index is inside this processor
972  if( (stress_system.solution->first_local_index() <= dof_index) &&
973  (dof_index < stress_system.solution->last_local_index()) )
974  {
975  stress_system.solution->set(dof_index, vonMises_value);
976  }
977  }
978 
979  // Close solution and update it
980  stress_system.solution->close();
981  stress_system.update();
982 };
set_homogeneous_physical_properties
void set_homogeneous_physical_properties(libMesh::EquationSystems &es, std::string &physicalParamsFile)
Set the homogeneous physical properties from a file.
Definition: assemble_functions_elasticity_3D.cpp:3
assemble_elasticity_with_weight_and_traction
void assemble_elasticity_with_weight_and_traction(libMesh::EquationSystems &es, const std::string &system_name, weight_parameter_function &weight_mask, WeightFunctionSystemType system_type, int traction_boundary_id, std::vector< double > traction_density)
Assemble homogeneous elasticity with domain weights and traction.
Definition: assemble_functions_elasticity_3D.cpp:397
weight_parameter_function
Definition: weight_parameter_function.h:13
eval_elasticity_tensor
libMesh::Real eval_elasticity_tensor(unsigned int i, unsigned int j, unsigned int k, unsigned int l, libMesh::Number E, libMesh::Number mu)
Calculate the elasticity tensor.
Definition: assemble_functions_elasticity_3D.cpp:180
set_heterogeneous_physical_properties
void set_heterogeneous_physical_properties(libMesh::EquationSystems &es, std::string &physicalParamsFile)
Set the heterogeneous, isotropic physical properties from a file.
Definition: assemble_functions_elasticity_3D.cpp:71
kronecker_delta
int kronecker_delta(unsigned int i, unsigned int j)
Definition: common_functions.cpp:3
compute_stresses
void compute_stresses(libMesh::EquationSystems &es)
Compute the stress (based on one of libMesh's examples)
Definition: assemble_functions_elasticity_3D.cpp:794
WeightFunctionSystemType
WeightFunctionSystemType
Enumerate used to define which weight function must be used to assemble the system.
Definition: ext_solver_libmesh_enums.h:14
assemble_elasticity_with_weight
void assemble_elasticity_with_weight(libMesh::EquationSystems &es, const std::string &system_name, weight_parameter_function &weight_mask, WeightFunctionSystemType system_type)
Assemble homogeneous elasticity with domain weights.
Definition: assemble_functions_elasticity_3D.cpp:223
assemble_elasticity_heterogeneous_with_weight
void assemble_elasticity_heterogeneous_with_weight(libMesh::EquationSystems &es, const std::string &system_name, weight_parameter_function &weight_mask, WeightFunctionSystemType system_type)
Assemble heterogeneous elasticity with domain weights.
Definition: assemble_functions_elasticity_3D.cpp:615
weight_parameter_function::get_alpha
double get_alpha(const libMesh::Point &qpoint, WeightFunctionSystemType system_type)
Definition: weight_parameter_function.h:137
eval_lambda_1
libMesh::Real eval_lambda_1(libMesh::Real E, libMesh::Real mu)
Calculate lambda_1 from E and Mu.
Definition: assemble_functions_elasticity_3D.cpp:175
assemble_functions_elasticity_3D.h
MASTER_bPerfLog_assemble_fem
const bool MASTER_bPerfLog_assemble_fem
Definition: assemble_functions_elasticity_3D.h:17
Update_SubK_isotropic
void Update_SubK_isotropic(libMesh::DenseSubMatrix< libMesh::Number > &SubK, unsigned int qp, unsigned int C_i, unsigned int C_k, const std::vector< std::vector< libMesh::RealGradient > > &dphi, const unsigned int n_components, const unsigned int n_u_dofs, const std::vector< libMesh::Real > &JxW, libMesh::Number E, libMesh::Number mu, double cte)
Calculate the rigidity sub-matrix contribution.
Definition: assemble_functions_elasticity_3D.cpp:195