constitutive_laws.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Non-inline functions for constitutive laws and strain-energy functions
27 
28 #include "constitutive_laws.h"
29 #include "../generic/elements.h"
30 
31 namespace oomph
32 {
33  //===============================================================
34  /// This function is used to check whether a matrix is square
35  //===============================================================
37  {
38  // If the number rows and columns is not equal, the matrix is not square
39  if (M.nrow() != M.ncol())
40  {
41  return false;
42  }
43  else
44  {
45  return true;
46  }
47  }
48 
49  //========================================================================
50  /// This function is used to check whether matrices are of equal dimension
51  //========================================================================
53  const DenseMatrix<double>& M1, const DenseMatrix<double>& M2)
54  {
55  // If the numbers of rows and columns are not the same, then the
56  // matrices are not of equal dimension
57  if ((M1.nrow() != M2.nrow()) || (M1.ncol() != M2.ncol()))
58  {
59  return false;
60  }
61  else
62  {
63  return true;
64  }
65  }
66 
67  //========================================================================
68  /// This function is used to provide simple error (bounce) checks on the
69  /// input to any calculate_second_piola_kirchhoff_stress
70  //=======================================================================
72  const DenseMatrix<double>& G,
73  DenseMatrix<double>& sigma)
74  {
75  // Test whether the undeformed metric tensor is square
76  if (!is_matrix_square(g))
77  {
78  throw OomphLibError("Undeformed metric tensor not square",
79  OOMPH_CURRENT_FUNCTION,
80  OOMPH_EXCEPTION_LOCATION);
81  }
82 
83  // If the deformed metric tensor does not have the same dimension as
84  // the undeformed tensor, complain
86  {
87  std::string error_message = "Deformed metric tensor does \n";
88  error_message +=
89  "not have same dimensions as the undeformed metric tensor\n";
90 
91  throw OomphLibError(
92  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
93  }
94 
95  // If the stress tensor does not have the same dimensions as the others
96  // complain.
97  if (!are_matrices_of_equal_dimensions(g, sigma))
98  {
99  std::string error_message =
100  "Strain tensor passed to calculate_green_strain() does \n";
101  error_message +=
102  "not have same dimensions as the undeformed metric tensor\n";
103 
104  throw OomphLibError(
105  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
106  }
107  }
108 
109 
110  //===========================================================================
111  /// The function to calculate the contravariant tensor from a covariant one
112  //===========================================================================
114  const DenseMatrix<double>& Gdown, DenseMatrix<double>& Gup)
115  {
116  // Initial error checking
117 #ifdef PARANOID
118  // Test that the matrices are of the same dimension
119  if (!are_matrices_of_equal_dimensions(Gdown, Gup))
120  {
121  throw OomphLibError("Matrices passed to calculate_contravariant() are "
122  "not of equal dimension",
123  OOMPH_CURRENT_FUNCTION,
124  OOMPH_EXCEPTION_LOCATION);
125  }
126 #endif
127 
128  // Find the dimension of the matrix
129  unsigned dim = Gdown.ncol();
130 
131  // If it's not square, I don't know what to do (yet)
132 #ifdef PARANOID
133  if (!is_matrix_square(Gdown))
134  {
135  std::string error_message =
136  "Matrix passed to calculate_contravariant() is not square\n";
137  error_message += "non-square matrix inversion not implemented yet!\n";
138 
139  throw OomphLibError(
140  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
141  }
142 #endif
143 
144  // Define the determinant of the matrix
145  double det = 0.0;
146 
147  // Now the inversion depends upon the dimension of the matrix
148  switch (dim)
149  {
150  // Zero dimensions
151  case 0:
152  throw OomphLibError("Zero dimensional matrix",
153  OOMPH_CURRENT_FUNCTION,
154  OOMPH_EXCEPTION_LOCATION);
155  break;
156 
157  // One dimension
158  case 1:
159  // The determinant is just the value of the only entry
160  det = Gdown(0, 0);
161  // The inverse is just the inverse of the value
162  Gup(0, 0) = 1.0 / Gdown(0, 0);
163  break;
164 
165 
166  // Two dimensions
167  case 2:
168  // Calculate the determinant
169  det = Gdown(0, 0) * Gdown(1, 1) - Gdown(0, 1) * Gdown(1, 0);
170  // Calculate entries of the contravariant tensor (inverse)
171  Gup(0, 0) = Gdown(1, 1) / det;
172  Gup(0, 1) = -Gdown(0, 1) / det;
173  Gup(1, 0) = -Gdown(1, 0) / det;
174  Gup(1, 1) = Gdown(0, 0) / det;
175  break;
176 
177  /// Three dimensions
178  case 3:
179  // Calculate the determinant of the matrix
180  det = Gdown(0, 0) * Gdown(1, 1) * Gdown(2, 2) +
181  Gdown(0, 1) * Gdown(1, 2) * Gdown(2, 0) +
182  Gdown(0, 2) * Gdown(1, 0) * Gdown(2, 1) -
183  Gdown(0, 0) * Gdown(1, 2) * Gdown(2, 1) -
184  Gdown(0, 1) * Gdown(1, 0) * Gdown(2, 2) -
185  Gdown(0, 2) * Gdown(1, 1) * Gdown(2, 0);
186 
187  // Calculate entries of the inverse matrix
188  Gup(0, 0) =
189  (Gdown(1, 1) * Gdown(2, 2) - Gdown(1, 2) * Gdown(2, 1)) / det;
190  Gup(0, 1) =
191  -(Gdown(0, 1) * Gdown(2, 2) - Gdown(0, 2) * Gdown(2, 1)) / det;
192  Gup(0, 2) =
193  (Gdown(0, 1) * Gdown(1, 2) - Gdown(0, 2) * Gdown(1, 1)) / det;
194  Gup(1, 0) =
195  -(Gdown(1, 0) * Gdown(2, 2) - Gdown(1, 2) * Gdown(2, 0)) / det;
196  Gup(1, 1) =
197  (Gdown(0, 0) * Gdown(2, 2) - Gdown(0, 2) * Gdown(2, 0)) / det;
198  Gup(1, 2) =
199  -(Gdown(0, 0) * Gdown(1, 2) - Gdown(0, 2) * Gdown(1, 0)) / det;
200  Gup(2, 0) =
201  (Gdown(1, 0) * Gdown(2, 1) - Gdown(1, 1) * Gdown(2, 0)) / det;
202  Gup(2, 1) =
203  -(Gdown(0, 0) * Gdown(2, 1) - Gdown(0, 1) * Gdown(2, 0)) / det;
204  Gup(2, 2) =
205  (Gdown(0, 0) * Gdown(1, 1) - Gdown(0, 1) * Gdown(1, 0)) / det;
206  break;
207 
208  default:
209  throw OomphLibError("Dimension of matrix must be 0, 1, 2 or 3\n",
210  OOMPH_CURRENT_FUNCTION,
211  OOMPH_EXCEPTION_LOCATION);
212  break;
213  }
214 
215  // Return the determinant of the matrix
216  return (det);
217  }
218 
219 
220  //===========================================================================
221  /// The function to calculate the derivatives of the contravariant tensor
222  /// and determinant of covariant tensor with respect to the components of
223  /// the covariant tensor
224  //===========================================================================
226  const DenseMatrix<double>& Gdown,
227  RankFourTensor<double>& d_Gup_dG,
228  DenseMatrix<double>& d_detG_dG)
229  {
230  // Find the dimension of the matrix
231  const unsigned dim = Gdown.ncol();
232 
233  // If it's not square, I don't know what to do (yet)
234 #ifdef PARANOID
235  if (!is_matrix_square(Gdown))
236  {
237  std::string error_message =
238  "Matrix passed to calculate_contravariant() is not square\n";
239  error_message += "non-square matrix inversion not implemented yet!\n";
240 
241  throw OomphLibError(
242  error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
243  }
244 #endif
245 
246  // Define the determinant of the matrix
247  double det = 0.0;
248 
249  // Now the inversion depends upon the dimension of the matrix
250  switch (dim)
251  {
252  // Zero dimensions
253  case 0:
254  throw OomphLibError("Zero dimensional matrix",
255  OOMPH_CURRENT_FUNCTION,
256  OOMPH_EXCEPTION_LOCATION);
257  break;
258 
259  // One dimension
260  case 1:
261  // There is only one entry, so derivatives are easy
262  d_detG_dG(0, 0) = 1.0;
263  d_Gup_dG(0, 0, 0, 0) = -1.0 / (Gdown(0, 0) * Gdown(0, 0));
264  break;
265 
266 
267  // Two dimensions
268  case 2:
269  // Calculate the determinant
270  det = Gdown(0, 0) * Gdown(1, 1) - Gdown(0, 1) * Gdown(1, 0);
271 
272  // Calculate the derivatives of the determinant
273  d_detG_dG(0, 0) = Gdown(1, 1);
274  // Need to use symmetry here
275  d_detG_dG(0, 1) = d_detG_dG(1, 0) = -2.0 * Gdown(0, 1);
276  d_detG_dG(1, 1) = Gdown(0, 0);
277 
278  // Calculate the "upper triangular" derivatives of the contravariant
279  // tensor
280  {
281  const double det2 = det * det;
282  d_Gup_dG(0, 0, 0, 0) = -Gdown(1, 1) * d_detG_dG(0, 0) / det2;
283  d_Gup_dG(0, 0, 0, 1) = -Gdown(1, 1) * d_detG_dG(0, 1) / det2;
284  d_Gup_dG(0, 0, 1, 1) =
285  1.0 / det - Gdown(1, 1) * d_detG_dG(1, 1) / det2;
286 
287  d_Gup_dG(0, 1, 0, 0) = Gdown(0, 1) * d_detG_dG(0, 0) / det2;
288  d_Gup_dG(0, 1, 0, 1) =
289  -1.0 / det + Gdown(0, 1) * d_detG_dG(0, 1) / det2;
290  d_Gup_dG(0, 1, 1, 1) = Gdown(0, 1) * d_detG_dG(1, 1) / det2;
291 
292  d_Gup_dG(1, 1, 0, 0) =
293  1.0 / det - Gdown(0, 0) * d_detG_dG(0, 0) / det2;
294  d_Gup_dG(1, 1, 0, 1) = -Gdown(0, 0) * d_detG_dG(0, 1) / det2;
295  d_Gup_dG(1, 1, 1, 1) = -Gdown(0, 0) * d_detG_dG(1, 1) / det2;
296  }
297 
298  // Calculate entries of the contravariant tensor (inverse)
299  // Gup(0,0) = Gdown(1,1)/det;
300  // Gup(0,1) = -Gdown(0,1)/det;
301  // Gup(1,0) = -Gdown(1,0)/det;
302  // Gup(1,1) = Gdown(0,0)/det;
303  break;
304 
305  /// Three dimensions
306  case 3:
307  // This is not yet implemented
308  throw OomphLibError(
309  "Analytic derivatives of metric tensors not yet implemented in 3D\n",
310  OOMPH_CURRENT_FUNCTION,
311  OOMPH_EXCEPTION_LOCATION);
312 
313  // Calculate the determinant of the matrix
314  det = Gdown(0, 0) * Gdown(1, 1) * Gdown(2, 2) +
315  Gdown(0, 1) * Gdown(1, 2) * Gdown(2, 0) +
316  Gdown(0, 2) * Gdown(1, 0) * Gdown(2, 1) -
317  Gdown(0, 0) * Gdown(1, 2) * Gdown(2, 1) -
318  Gdown(0, 1) * Gdown(1, 0) * Gdown(2, 2) -
319  Gdown(0, 2) * Gdown(1, 1) * Gdown(2, 0);
320 
321  // Calculate entries of the inverse matrix
322  // Gup(0,0) = (Gdown(1,1)*Gdown(2,2) - Gdown(1,2)*Gdown(2,1))/det;
323  // Gup(0,1) = -(Gdown(0,1)*Gdown(2,2) - Gdown(0,2)*Gdown(2,1))/det;
324  // Gup(0,2) = (Gdown(0,1)*Gdown(1,2) - Gdown(0,2)*Gdown(1,1))/det;
325  // Gup(1,0) = -(Gdown(1,0)*Gdown(2,2) - Gdown(1,2)*Gdown(2,0))/det;
326  // Gup(1,1) = (Gdown(0,0)*Gdown(2,2) - Gdown(0,2)*Gdown(2,0))/det;
327  // Gup(1,2) = -(Gdown(0,0)*Gdown(1,2) - Gdown(0,2)*Gdown(1,0))/det;
328  // Gup(2,0) = (Gdown(1,0)*Gdown(2,1) - Gdown(1,1)*Gdown(2,0))/det;
329  // Gup(2,1) = -(Gdown(0,0)*Gdown(2,1) - Gdown(0,1)*Gdown(2,0))/det;
330  // Gup(2,2) = (Gdown(0,0)*Gdown(1,1) - Gdown(0,1)*Gdown(1,0))/det;
331  break;
332 
333  default:
334  throw OomphLibError("Dimension of matrix must be 0, 1, 2 or 3\n",
335  OOMPH_CURRENT_FUNCTION,
336  OOMPH_EXCEPTION_LOCATION);
337  break;
338  }
339  }
340 
341 
342  //=========================================================================
343  /// Calculate the derivatives of the contravariant
344  /// 2nd Piola Kirchhoff stress tensor with respect to the deformed metric
345  /// tensor. Arguments are the
346  /// covariant undeformed and deformed metric tensor and the
347  /// matrix in which to return the derivatives of the stress tensor
348  /// The default implementation uses finite differences, but can be
349  /// overloaded for constitutive laws in which an analytic formulation
350  /// is possible.
351  //==========================================================================
353  const DenseMatrix<double>& g,
354  const DenseMatrix<double>& G,
355  const DenseMatrix<double>& sigma,
356  RankFourTensor<double>& d_sigma_dG,
357  const bool& symmetrize_tensor)
358  {
359  // Initial error checking
360 #ifdef PARANOID
361  // Test that the matrices are of the same dimension
363  {
364  throw OomphLibError("Matrices passed are not of equal dimension",
365  OOMPH_CURRENT_FUNCTION,
366  OOMPH_EXCEPTION_LOCATION);
367  }
368 #endif
369 
370  // Find the dimension of the matrix (assuming that it's square)
371  const unsigned dim = G.ncol();
372 
373  // Find the dimension
374  // FD step
375  const double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
376 
377  // Advanced metric tensor
378  DenseMatrix<double> G_pls(dim, dim);
379  DenseMatrix<double> sigma_pls(dim, dim);
380 
381  // Copy across the original value
382  for (unsigned i = 0; i < dim; i++)
383  {
384  for (unsigned j = 0; j < dim; j++)
385  {
386  G_pls(i, j) = G(i, j);
387  }
388  }
389 
390  // Do FD -- only w.r.t. to upper indices, exploiting symmetry.
391  // NOTE: We exploit the symmetry of the stress and metric tensors
392  // by incrementing G(i,j) and G(j,i) simultaenously and
393  // only fill in the "upper" triangles without copying things
394  // across the lower triangle. This is taken into account
395  // in the solid mechanics codes.
396  for (unsigned i = 0; i < dim; i++)
397  {
398  for (unsigned j = i; j < dim; j++)
399  {
400  G_pls(i, j) += eps_fd;
401  G_pls(j, i) = G_pls(i, j);
402 
403  // Get advanced stress
404  this->calculate_second_piola_kirchhoff_stress(g, G_pls, sigma_pls);
405 
406  for (unsigned ii = 0; ii < dim; ii++)
407  {
408  for (unsigned jj = ii; jj < dim; jj++)
409  {
410  d_sigma_dG(ii, jj, i, j) =
411  (sigma_pls(ii, jj) - sigma(ii, jj)) / eps_fd;
412  }
413  }
414 
415  // Reset
416  G_pls(i, j) = G(i, j);
417  G_pls(j, i) = G(j, i);
418  }
419  }
420 
421  // If we are symmetrising the tensor, do so
422  if (symmetrize_tensor)
423  {
424  for (unsigned i = 0; i < dim; i++)
425  {
426  for (unsigned j = 0; j < i; j++)
427  {
428  for (unsigned ii = 0; ii < dim; ii++)
429  {
430  for (unsigned jj = 0; jj < ii; jj++)
431  {
432  d_sigma_dG(ii, jj, i, j) = d_sigma_dG(jj, ii, j, i);
433  }
434  }
435  }
436  }
437  }
438  }
439 
440 
441  //=========================================================================
442  /// Calculate the derivatives of the contravariant
443  /// 2nd Piola Kirchhoff stress tensor \f$ \sigma^{ij}\f$.
444  /// with respect to the deformed metric tensor.
445  /// Also return the derivatives of the determinant of the
446  /// deformed metric tensor with respect to the deformed metric tensor.
447  /// This form is appropriate
448  /// for truly-incompressible materials.
449  /// The default implementation uses finite differences for the
450  /// derivatives that depend on the constitutive law, but not
451  /// for the derivatives of the determinant, which are generic.
452  //========================================================================
454  const DenseMatrix<double>& g,
455  const DenseMatrix<double>& G,
456  const DenseMatrix<double>& sigma,
457  const double& detG,
458  const double& interpolated_solid_p,
459  RankFourTensor<double>& d_sigma_dG,
460  DenseMatrix<double>& d_detG_dG,
461  const bool& symmetrize_tensor)
462  {
463  // Initial error checking
464 #ifdef PARANOID
465  // Test that the matrices are of the same dimension
467  {
468  throw OomphLibError("Matrices passed are not of equal dimension",
469  OOMPH_CURRENT_FUNCTION,
470  OOMPH_EXCEPTION_LOCATION);
471  }
472 #endif
473 
474  // Find the dimension of the matrix (assuming that it's square)
475  const unsigned dim = G.ncol();
476 
477  // FD step
478  const double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
479 
480  // Advanced metric tensor etc.
481  DenseMatrix<double> G_pls(dim, dim);
482  DenseMatrix<double> sigma_dev_pls(dim, dim);
483  DenseMatrix<double> Gup_pls(dim, dim);
484  double detG_pls;
485 
486  // Copy across
487  for (unsigned i = 0; i < dim; i++)
488  {
489  for (unsigned j = 0; j < dim; j++)
490  {
491  G_pls(i, j) = G(i, j);
492  }
493  }
494 
495  // Do FD -- only w.r.t. to upper indices, exploiting symmetry.
496  // NOTE: We exploit the symmetry of the stress and metric tensors
497  // by incrementing G(i,j) and G(j,i) simultaenously and
498  // only fill in the "upper" triangles without copying things
499  // across the lower triangle. This is taken into account
500  // in the remaining code further below.
501  for (unsigned i = 0; i < dim; i++)
502  {
503  for (unsigned j = i; j < dim; j++)
504  {
505  G_pls(i, j) += eps_fd;
506  G_pls(j, i) = G_pls(i, j);
507 
508  // Get advanced stress
510  g, G_pls, sigma_dev_pls, Gup_pls, detG_pls);
511 
512 
513  // Derivative of determinant of deformed metric tensor
514  d_detG_dG(i, j) = (detG_pls - detG) / eps_fd;
515 
516  // Derivatives of deviatoric stress and "upper" deformed metric
517  // tensor
518  for (unsigned ii = 0; ii < dim; ii++)
519  {
520  for (unsigned jj = ii; jj < dim; jj++)
521  {
522  d_sigma_dG(ii, jj, i, j) =
523  (sigma_dev_pls(ii, jj) - interpolated_solid_p * Gup_pls(ii, jj) -
524  sigma(ii, jj)) /
525  eps_fd;
526  }
527  }
528 
529  // Reset
530  G_pls(i, j) = G(i, j);
531  G_pls(j, i) = G(j, i);
532  }
533  }
534 
535  // If we are symmetrising the tensor, do so
536  if (symmetrize_tensor)
537  {
538  for (unsigned i = 0; i < dim; i++)
539  {
540  for (unsigned j = 0; j < i; j++)
541  {
542  d_detG_dG(i, j) = d_detG_dG(j, i);
543 
544  for (unsigned ii = 0; ii < dim; ii++)
545  {
546  for (unsigned jj = 0; jj < ii; jj++)
547  {
548  d_sigma_dG(ii, jj, i, j) = d_sigma_dG(jj, ii, j, i);
549  }
550  }
551  }
552  }
553  }
554  }
555 
556 
557  //========================================================================
558  /// Calculate the derivatives of the contravariant
559  /// 2nd Piola Kirchoff stress tensor with respect to the deformed metric
560  /// tensor. Also return the derivatives of the generalised dilatation,
561  /// \f$ d, \f$ with respect to the deformed metric tensor.
562  /// This form is appropriate
563  /// for near-incompressible materials.
564  /// The default implementation uses finite differences.
565  //=======================================================================
567  const DenseMatrix<double>& g,
568  const DenseMatrix<double>& G,
569  const DenseMatrix<double>& sigma,
570  const double& gen_dil,
571  const double& inv_kappa,
572  const double& interpolated_solid_p,
573  RankFourTensor<double>& d_sigma_dG,
574  DenseMatrix<double>& d_gen_dil_dG,
575  const bool& symmetrize_tensor)
576  {
577  // Initial error checking
578 #ifdef PARANOID
579  // Test that the matrices are of the same dimension
581  {
582  throw OomphLibError("Matrices passed are not of equal dimension",
583  OOMPH_CURRENT_FUNCTION,
584  OOMPH_EXCEPTION_LOCATION);
585  }
586 #endif
587 
588  // Find the dimension of the matrix (assuming that it's square)
589  const unsigned dim = G.ncol();
590 
591  // FD step
592  const double eps_fd = GeneralisedElement::Default_fd_jacobian_step;
593 
594  // Advanced metric tensor etc
595  DenseMatrix<double> G_pls(dim, dim);
596  DenseMatrix<double> sigma_dev_pls(dim, dim);
597  DenseMatrix<double> Gup_pls(dim, dim);
598  double gen_dil_pls;
599  double inv_kappa_pls;
600 
601  // Copy across
602  for (unsigned i = 0; i < dim; i++)
603  {
604  for (unsigned j = 0; j < dim; j++)
605  {
606  G_pls(i, j) = G(i, j);
607  }
608  }
609 
610  // Do FD -- only w.r.t. to upper indices, exploiting symmetry.
611  // NOTE: We exploit the symmetry of the stress and metric tensors
612  // by incrementing G(i,j) and G(j,i) simultaenously and
613  // only fill in the "upper" triangles without copying things
614  // across the lower triangle. This is taken into account
615  // in the remaining code further below.
616  for (unsigned i = 0; i < dim; i++)
617  {
618  for (unsigned j = i; j < dim; j++)
619  {
620  G_pls(i, j) += eps_fd;
621  G_pls(j, i) = G_pls(i, j);
622 
623  // Get advanced stress
625  g, G_pls, sigma_dev_pls, Gup_pls, gen_dil_pls, inv_kappa_pls);
626 
627  // Derivative of generalised dilatation
628  d_gen_dil_dG(i, j) = (gen_dil_pls - gen_dil) / eps_fd;
629 
630  // Derivatives of deviatoric stress and "upper" deformed metric
631  // tensor
632  for (unsigned ii = 0; ii < dim; ii++)
633  {
634  for (unsigned jj = ii; jj < dim; jj++)
635  {
636  d_sigma_dG(ii, jj, i, j) =
637  (sigma_dev_pls(ii, jj) - interpolated_solid_p * Gup_pls(ii, jj) -
638  sigma(ii, jj)) /
639  eps_fd;
640  }
641  }
642 
643  // Reset
644  G_pls(i, j) = G(i, j);
645  G_pls(j, i) = G(j, i);
646  }
647  }
648 
649  // If we are symmetrising the tensor, do so
650  if (symmetrize_tensor)
651  {
652  for (unsigned i = 0; i < dim; i++)
653  {
654  for (unsigned j = 0; j < i; j++)
655  {
656  d_gen_dil_dG(i, j) = d_gen_dil_dG(j, i);
657 
658  for (unsigned ii = 0; ii < dim; ii++)
659  {
660  for (unsigned jj = 0; jj < ii; jj++)
661  {
662  d_sigma_dG(ii, jj, i, j) = d_sigma_dG(jj, ii, j, i);
663  }
664  }
665  }
666  }
667  }
668  }
669 
670 
671  /// //////////////////////////////////////////////////////////////////
672  /// //////////////////////////////////////////////////////////////////
673  /// //////////////////////////////////////////////////////////////////
674 
675 
676  //=====================================================================
677  /// Calculate the contravariant 2nd Piola Kirchhoff
678  /// stress tensor. Arguments are the
679  /// covariant undeformed (stress-free) and deformed metric
680  /// tensors, g and G, and the matrix in which to return the stress tensor.
681  //=====================================================================
683  const DenseMatrix<double>& g,
684  const DenseMatrix<double>& G,
685  DenseMatrix<double>& sigma)
686  {
687  // Error checking
688 #ifdef PARANOID
689  error_checking_in_input(g, G, sigma);
690 #endif
691 
692  // Find the dimension of the problem
693  unsigned dim = G.nrow();
694 
695  // Calculate the contravariant Deformed metric tensor
696  DenseMatrix<double> Gup(dim);
697  // We don't need the Jacobian so cast the function to void
698  (void)calculate_contravariant(G, Gup);
699 
700  // Premultiply some constants
701  double C1 = (*E_pt) / (2.0 * (1.0 + (*Nu_pt)));
702  double C2 = 2.0 * (*Nu_pt) / (1.0 - 2.0 * (*Nu_pt));
703 
704  // Strain tensor
705  DenseMatrix<double> strain(dim, dim);
706 
707  // Upper triangle
708  for (unsigned i = 0; i < dim; i++)
709  {
710  for (unsigned j = i; j < dim; j++)
711  {
712  strain(i, j) = 0.5 * (G(i, j) - g(i, j));
713  }
714  }
715 
716  // Copy across
717  for (unsigned i = 0; i < dim; i++)
718  {
719  for (unsigned j = 0; j < i; j++)
720  {
721  strain(i, j) = strain(j, i);
722  }
723  }
724 
725  // Compute upper triangle of stress
726  for (unsigned i = 0; i < dim; i++)
727  {
728  for (unsigned j = i; j < dim; j++)
729  {
730  // Initialise this component of sigma
731  sigma(i, j) = 0.0;
732  for (unsigned k = 0; k < dim; k++)
733  {
734  for (unsigned l = 0; l < dim; l++)
735  {
736  sigma(i, j) += C1 *
737  (Gup(i, k) * Gup(j, l) + Gup(i, l) * Gup(j, k) +
738  C2 * Gup(i, j) * Gup(k, l)) *
739  strain(k, l);
740  }
741  }
742  }
743  }
744 
745  // Copy across
746  for (unsigned i = 0; i < dim; i++)
747  {
748  for (unsigned j = 0; j < i; j++)
749  {
750  sigma(i, j) = sigma(j, i);
751  }
752  }
753  }
754 
755  //===========================================================================
756  /// Calculate the deviatoric part of the contravariant
757  /// 2nd Piola Kirchoff stress tensor. Also return the contravariant
758  /// deformed metric tensor, the generalised dilatation, \f$ d, \f$ and
759  /// the inverse of the bulk modulus \f$ \kappa\f$. This form is appropriate
760  /// for near-incompressible materials for which
761  /// \f$ \sigma^{ij} = -p G^{ij} + \overline{ \sigma^{ij}} \f$
762  /// where the "pressure" \f$ p \f$ is determined from
763  /// \f$ p / \kappa - d =0 \f$.
764  //===========================================================================
766  const DenseMatrix<double>& g,
767  const DenseMatrix<double>& G,
768  DenseMatrix<double>& sigma_dev,
769  DenseMatrix<double>& Gup,
770  double& gen_dil,
771  double& inv_kappa)
772  {
773  // Find the dimension of the problem
774  unsigned dim = G.nrow();
775 
776  // Assign memory for the determinant of the deformed metric tensor
777  double detG = 0.0;
778 
779  // Compute deviatoric stress by calling the incompressible
780  // version of this function
781  calculate_second_piola_kirchhoff_stress(g, G, sigma_dev, Gup, detG);
782 
783  // Calculate the inverse of the "bulk" modulus
784  inv_kappa =
785  (1.0 - 2.0 * (*Nu_pt)) * (1.0 + (*Nu_pt)) / ((*E_pt) * (*Nu_pt));
786 
787  // Finally compute the generalised dilatation (i.e. the term that
788  // must be zero if \kappa \to \infty
789  gen_dil = 0.0;
790  for (unsigned i = 0; i < dim; i++)
791  {
792  for (unsigned j = 0; j < dim; j++)
793  {
794  gen_dil += Gup(i, j) * 0.5 * (G(i, j) - g(i, j));
795  }
796  }
797  }
798 
799  //======================================================================
800  /// Calculate the deviatoric part
801  /// \f$ \overline{ \sigma^{ij}}\f$ of the contravariant
802  /// 2nd Piola Kirchhoff stress tensor \f$ \sigma^{ij}\f$.
803  /// Also return the contravariant deformed metric
804  /// tensor and the determinant of the deformed metric tensor.
805  /// This form is appropriate
806  /// for truly-incompressible materials for which
807  /// \f$ \sigma^{ij} = - p G^{ij} +\overline{ \sigma^{ij}} \f$
808  /// where the "pressure" \f$ p \f$ is determined by
809  /// \f$ \det G_{ij} - \det g_{ij} = 0 \f$.
810  //======================================================================
812  const DenseMatrix<double>& g,
813  const DenseMatrix<double>& G,
814  DenseMatrix<double>& sigma_dev,
815  DenseMatrix<double>& Gup,
816  double& detG)
817  {
818  // Error checking
819 #ifdef PARANOID
820  error_checking_in_input(g, G, sigma_dev);
821 #endif
822 
823  // Find the dimension of the problem
824  unsigned dim = G.nrow();
825 
826  // Calculate the contravariant Deformed metric tensor
827  detG = calculate_contravariant(G, Gup);
828 
829  // Premultiply the appropriate physical constant
830  double C1 = (*E_pt) / (2.0 * (1.0 + (*Nu_pt)));
831 
832  // Strain tensor
833  DenseMatrix<double> strain(dim, dim);
834 
835  // Upper triangle
836  for (unsigned i = 0; i < dim; i++)
837  {
838  for (unsigned j = i; j < dim; j++)
839  {
840  strain(i, j) = 0.5 * (G(i, j) - g(i, j));
841  }
842  }
843 
844  // Copy across
845  for (unsigned i = 0; i < dim; i++)
846  {
847  for (unsigned j = 0; j < i; j++)
848  {
849  strain(i, j) = strain(j, i);
850  }
851  }
852 
853  // Compute upper triangle of stress
854  for (unsigned i = 0; i < dim; i++)
855  {
856  for (unsigned j = i; j < dim; j++)
857  {
858  // Initialise this component of sigma
859  sigma_dev(i, j) = 0.0;
860  for (unsigned k = 0; k < dim; k++)
861  {
862  for (unsigned l = 0; l < dim; l++)
863  {
864  sigma_dev(i, j) += C1 *
865  (Gup(i, k) * Gup(j, l) + Gup(i, l) * Gup(j, k)) *
866  strain(k, l);
867  }
868  }
869  }
870  }
871 
872  // Copy across
873  for (unsigned i = 0; i < dim; i++)
874  {
875  for (unsigned j = 0; j < i; j++)
876  {
877  sigma_dev(i, j) = sigma_dev(j, i);
878  }
879  }
880  }
881 
882 
883  /// //////////////////////////////////////////////////////////////////////
884  /// //////////////////////////////////////////////////////////////////////
885  /// //////////////////////////////////////////////////////////////////////
886 
887 
888  //========================================================================
889  /// Calculate the contravariant 2nd Piola Kirchhoff
890  /// stress tensor. Arguments are the
891  /// covariant undeformed and deformed metric tensor and the
892  /// matrix in which to return the stress tensor.
893  /// Uses correct 3D invariants for 2D (plane strain) problems.
894  //=======================================================================
897  const DenseMatrix<double>& G,
898  DenseMatrix<double>& sigma)
899  {
900 // Error checking
901 #ifdef PARANOID
902  error_checking_in_input(g, G, sigma);
903 #endif
904 
905  // Find the dimension of the problem
906  unsigned dim = g.nrow();
907 
908 #ifdef PARANOID
909  if (dim == 1)
910  {
911  std::string function_name =
912  "IsotropicStrainEnergyFunctionConstitutiveLaw::";
913  function_name += "calculate_second_piola_kirchhoff_stress()";
914 
915  throw OomphLibError("Check constitutive equations carefully when dim=1",
916  OOMPH_CURRENT_FUNCTION,
917  OOMPH_EXCEPTION_LOCATION);
918  }
919 #endif
920 
921  // Calculate the contravariant undeformed and deformed metric tensors
922  // and get the determinants of the metric tensors
923  DenseMatrix<double> gup(dim), Gup(dim);
924  double detg = calculate_contravariant(g, gup);
925  double detG = calculate_contravariant(G, Gup);
926 
927  // Calculate the strain invariants
928  Vector<double> I(3, 0.0);
929  // The third strain invaraint is the volumetric change
930  I[2] = detG / detg;
931  // The first and second are a bit more complex --- see G&Z
932  for (unsigned i = 0; i < dim; i++)
933  {
934  for (unsigned j = 0; j < dim; j++)
935  {
936  I[0] += gup(i, j) * G(i, j);
937  I[1] += g(i, j) * Gup(i, j);
938  }
939  }
940 
941  // If 2D we assume plane strain: In this case the 3D tensors have
942  // a 1 on the diagonal and zeroes in the off-diagonals of their
943  // third rows and columns. Only effect: Increase the first two
944  // invariants by one; rest of the computation can just be performed
945  // over the 2d set of coordinates.
946  if (dim == 2)
947  {
948  I[0] += 1.0;
949  I[1] += 1.0;
950  }
951 
952  // Second strain invariant is multiplied by the third.
953  I[1] *= I[2];
954 
955  // Calculate the derivatives of the strain energy function wrt the
956  // strain invariants
957  Vector<double> dWdI(3, 0.0);
959 
960 
961  // Only bother to compute the tensor B^{ij} (Green & Zerna notation)
962  // if the derivative wrt the second strain invariant is non-zero
963  DenseMatrix<double> Bup(dim, dim, 0.0);
964  if (std::fabs(dWdI[1]) > 0.0)
965  {
966  for (unsigned i = 0; i < dim; i++)
967  {
968  for (unsigned j = 0; j < dim; j++)
969  {
970  Bup(i, j) = I[0] * gup(i, j);
971  for (unsigned r = 0; r < dim; r++)
972  {
973  for (unsigned s = 0; s < dim; s++)
974  {
975  Bup(i, j) -= gup(i, r) * gup(j, s) * G(r, s);
976  }
977  }
978  }
979  }
980  }
981 
982  // Now set the values of the functions phi, psi and p (Green & Zerna
983  // notation) Note that the Green & Zerna stress \tau^{ij} is
984  // s^{ij}/sqrt(I[2]), where s^{ij} is the desired second Piola-Kirchhoff
985  // stress tensor so we multiply their constants by sqrt(I[2])
986  double phi = 2.0 * dWdI[0];
987  double psi = 2.0 * dWdI[1];
988  double p = 2.0 * dWdI[2] * I[2];
989 
990  // Put it all together to get the stress
991  for (unsigned i = 0; i < dim; i++)
992  {
993  for (unsigned j = 0; j < dim; j++)
994  {
995  sigma(i, j) = phi * gup(i, j) + psi * Bup(i, j) + p * Gup(i, j);
996  }
997  }
998  }
999 
1000  //===========================================================================
1001  /// Calculate the deviatoric part
1002  /// \f$ \overline{ \sigma^{ij}}\f$ of the contravariant
1003  /// 2nd Piola Kirchhoff stress tensor \f$ \sigma^{ij}\f$.
1004  /// Also return the contravariant deformed metric
1005  /// tensor and the determinant of the deformed metric tensor.
1006  /// Uses correct 3D invariants for 2D (plane strain) problems.
1007  /// This is the version for the pure incompressible formulation.
1008  //============================================================================
1011  const DenseMatrix<double>& G,
1012  DenseMatrix<double>& sigma_dev,
1013  DenseMatrix<double>& Gup,
1014  double& detG)
1015  {
1016 // Error checking
1017 #ifdef PARANOID
1018  error_checking_in_input(g, G, sigma_dev);
1019 #endif
1020 
1021  // Find the dimension of the problem
1022  unsigned dim = g.nrow();
1023 
1024 
1025 #ifdef PARANOID
1026  if (dim == 1)
1027  {
1028  std::string function_name =
1029  "IsotropicStrainEnergyFunctionConstitutiveLaw::";
1030  function_name += "calculate_second_piola_kirchhoff_stress()";
1031 
1032  throw OomphLibError("Check constitutive equations carefully when dim=1",
1033  OOMPH_CURRENT_FUNCTION,
1034  OOMPH_EXCEPTION_LOCATION);
1035  }
1036 #endif
1037 
1038  // Calculate the contravariant undeformed and deformed metric tensors
1039  DenseMatrix<double> gup(dim);
1040  // Don't need this determinant
1041  (void)calculate_contravariant(g, gup);
1042  // These are passed back
1043  detG = calculate_contravariant(G, Gup);
1044 
1045  // Calculate the strain invariants
1046  Vector<double> I(3, 0.0);
1047  // The third strain invaraint must be one (incompressibility)
1048  I[2] = 1.0;
1049  // The first and second are a bit more complex
1050  for (unsigned i = 0; i < dim; i++)
1051  {
1052  for (unsigned j = 0; j < dim; j++)
1053  {
1054  I[0] += gup(i, j) * G(i, j);
1055  I[1] += g(i, j) * Gup(i, j);
1056  }
1057  }
1058 
1059  // If 2D we assume plane strain: In this case the 3D tensors have
1060  // a 1 on the diagonal and zeroes in the off-diagonals of their
1061  // third rows and columns. Only effect: Increase the first two
1062  // invariants by one; rest of the computation can just be performed
1063  // over the 2d set of coordinates.
1064  if (dim == 2)
1065  {
1066  I[0] += 1.0;
1067  I[1] += 1.0;
1068  }
1069 
1070  // Calculate the derivatives of the strain energy function wrt the
1071  // strain invariants
1072  Vector<double> dWdI(3, 0.0);
1074 
1075  // Only bother to compute the tensor B^{ij} (Green & Zerna notation)
1076  // if the derivative wrt the second strain invariant is non-zero
1077  DenseMatrix<double> Bup(dim, dim, 0.0);
1078  if (std::fabs(dWdI[1]) > 0.0)
1079  {
1080  for (unsigned i = 0; i < dim; i++)
1081  {
1082  for (unsigned j = 0; j < dim; j++)
1083  {
1084  Bup(i, j) = I[0] * gup(i, j);
1085  for (unsigned r = 0; r < dim; r++)
1086  {
1087  for (unsigned s = 0; s < dim; s++)
1088  {
1089  Bup(i, j) -= gup(i, r) * gup(j, s) * G(r, s);
1090  }
1091  }
1092  }
1093  }
1094  }
1095 
1096  // Now set the values of the functions phi and psi (Green & Zerna notation)
1097  double phi = 2.0 * dWdI[0];
1098  double psi = 2.0 * dWdI[1];
1099  // Calculate the trace/dim of the first two terms of the stress tensor
1100  // phi g^{ij} + psi B^{ij} (see Green & Zerna)
1101  double K;
1102  // In two-d, we cannot use the strain invariants directly
1103  // but can use symmetry of the tensors involved
1104  if (dim == 2)
1105  {
1106  K = 0.5 * ((I[0] - 1.0) * phi +
1107  psi * (Bup(0, 0) * G(0, 0) + Bup(1, 1) * G(1, 1) +
1108  2.0 * Bup(0, 1) * G(0, 1)));
1109  }
1110  // In three-d we can make use of the strain invariants, see Green & Zerna
1111  else
1112  {
1113  K = (I[0] * phi + 2.0 * I[1] * psi) / 3.0;
1114  }
1115 
1116  // Put it all together to get the stress, subtracting the trace of the
1117  // first two terms to ensure that the stress is deviatoric, which means
1118  // that the computed pressure is the mechanical pressure
1119  for (unsigned i = 0; i < dim; i++)
1120  {
1121  for (unsigned j = 0; j < dim; j++)
1122  {
1123  sigma_dev(i, j) = phi * gup(i, j) + psi * Bup(i, j) - K * Gup(i, j);
1124  }
1125  }
1126  }
1127 
1128  //===========================================================================
1129  /// Calculate the deviatoric part of the contravariant
1130  /// 2nd Piola Kirchoff stress tensor. Also return the contravariant
1131  /// deformed metric tensor, the generalised dilatation, \f$ d, \f$ and
1132  /// the inverse of the bulk modulus \f$ \kappa\f$.
1133  /// Uses correct 3D invariants for 2D (plane strain) problems.
1134  /// This is the version for the near-incompressible formulation.
1135  //===========================================================================
1138  const DenseMatrix<double>& G,
1139  DenseMatrix<double>& sigma_dev,
1140  DenseMatrix<double>& Gup,
1141  double& gen_dil,
1142  double& inv_kappa)
1143  {
1144 // Error checking
1145 #ifdef PARANOID
1146  error_checking_in_input(g, G, sigma_dev);
1147 #endif
1148 
1149  // Find the dimension of the problem
1150  unsigned dim = g.nrow();
1151 
1152 #ifdef PARANOID
1153  if (dim == 1)
1154  {
1155  std::string function_name =
1156  "IsotropicStrainEnergyFunctionConstitutiveLaw::";
1157  function_name += "calculate_second_piola_kirchhoff_stress()";
1158 
1159  throw OomphLibError("Check constitutive equations carefully when dim=1",
1160  OOMPH_CURRENT_FUNCTION,
1161  OOMPH_EXCEPTION_LOCATION);
1162  }
1163 #endif
1164 
1165  // Calculate the contravariant undeformed and deformed metric tensors
1166  // and get the determinants of the metric tensors
1167  DenseMatrix<double> gup(dim);
1168  double detg = calculate_contravariant(g, gup);
1169  double detG = calculate_contravariant(G, Gup);
1170 
1171  // Calculate the strain invariants
1172  Vector<double> I(3, 0.0);
1173  // The third strain invaraint is the volumetric change
1174  I[2] = detG / detg;
1175  // The first and second are a bit more complex --- see G&Z
1176  for (unsigned i = 0; i < dim; i++)
1177  {
1178  for (unsigned j = 0; j < dim; j++)
1179  {
1180  I[0] += gup(i, j) * G(i, j);
1181  I[1] += g(i, j) * Gup(i, j);
1182  }
1183  }
1184 
1185  // If 2D we assume plane strain: In this case the 3D tensors have
1186  // a 1 on the diagonal and zeroes in the off-diagonals of their
1187  // third rows and columns. Only effect: Increase the first two
1188  // invariants by one; rest of the computation can just be performed
1189  // over the 2d set of coordinates.
1190  if (dim == 2)
1191  {
1192  I[0] += 1.0;
1193  I[1] += 1.0;
1194  }
1195 
1196  // Second strain invariant is multiplied by the third.
1197  I[1] *= I[2];
1198 
1199  // Calculate the derivatives of the strain energy function wrt the
1200  // strain invariants
1201  Vector<double> dWdI(3, 0.0);
1203 
1204  // Only bother to calculate the tensor B^{ij} (Green & Zerna notation)
1205  // if the derivative wrt the second strain invariant is non-zero
1206  DenseMatrix<double> Bup(dim, dim, 0.0);
1207  if (std::fabs(dWdI[1]) > 0.0)
1208  {
1209  for (unsigned i = 0; i < dim; i++)
1210  {
1211  for (unsigned j = 0; j < dim; j++)
1212  {
1213  Bup(i, j) = I[0] * gup(i, j);
1214  for (unsigned r = 0; r < dim; r++)
1215  {
1216  for (unsigned s = 0; s < dim; s++)
1217  {
1218  Bup(i, j) -= gup(i, r) * gup(j, s) * G(r, s);
1219  }
1220  }
1221  }
1222  }
1223  }
1224 
1225  // Now set the values of the functions phi and psi (Green & Zerna notation)
1226  // but multiplied by sqrt(I[2]) to recover the second Piola-Kirchhoff stress
1227  double phi = 2.0 * dWdI[0];
1228  double psi = 2.0 * dWdI[1];
1229 
1230  // Calculate the trace/dim of the first two terms of the stress tensor
1231  // phi g^{ij} + psi B^{ij} (see Green & Zerna)
1232  double K;
1233  // In two-d, we cannot use the strain invariants directly,
1234  // but we can use symmetry of the tensors involved
1235  if (dim == 2)
1236  {
1237  K = 0.5 * ((I[0] - 1.0) * phi +
1238  psi * (Bup(0, 0) * G(0, 0) + Bup(1, 1) * G(1, 1) +
1239  2.0 * Bup(0, 1) * G(0, 1)));
1240  }
1241  // In three-d we can make use of the strain invariants
1242  else
1243  {
1244  K = (I[0] * phi + 2.0 * I[1] * psi) / 3.0;
1245  }
1246 
1247  // Choose inverse kappa to be one...
1248  inv_kappa = 1.0;
1249 
1250  //...then the generalised dilation is the same as p in Green & Zerna's
1251  // notation, but multiplied by sqrt(I[2]) with the addition of the
1252  // terms that are subtracted to make the other part of the stress
1253  // deviatoric
1254  gen_dil = 2.0 * dWdI[2] * I[2] + K;
1255 
1256  // Calculate the deviatoric part of the stress by subtracting
1257  // the computed trace/dim
1258  for (unsigned i = 0; i < dim; i++)
1259  {
1260  for (unsigned j = 0; j < dim; j++)
1261  {
1262  sigma_dev(i, j) = phi * gup(i, j) + psi * Bup(i, j) - K * Gup(i, j);
1263  }
1264  }
1265  }
1266 
1267 } // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
void error_checking_in_input(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
Check for errors in the input, i.e. check that the dimensions of the arrays are all consistent.
double calculate_contravariant(const DenseMatrix< double > &Gcov, DenseMatrix< double > &Gcontra)
Calculate a contravariant tensor from a covariant tensor, and return the determinant of the covariant...
virtual void calculate_d_second_piola_kirchhoff_stress_dG(const DenseMatrix< double > &g, const DenseMatrix< double > &G, const DenseMatrix< double > &sigma, RankFourTensor< double > &d_sigma_dG, const bool &symmetrize_tensor=true)
Calculate the derivatives of the contravariant 2nd Piola Kirchhoff stress tensor with respect to the ...
void calculate_d_contravariant_dG(const DenseMatrix< double > &Gcov, RankFourTensor< double > &dGcontra_dG, DenseMatrix< double > &d_detG_dG)
Calculate the derivatives of the contravariant tensor and the derivatives of the determinant of the c...
bool is_matrix_square(const DenseMatrix< double > &M)
Test whether a matrix is square.
bool are_matrices_of_equal_dimensions(const DenseMatrix< double > &M1, const DenseMatrix< double > &M2)
Test whether two matrices are of equal dimensions.
virtual void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)=0
Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed...
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:485
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:491
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition: elements.h:1198
void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed...
double * E_pt
Young's modulus.
double * Nu_pt
Poisson ratio.
StrainEnergyFunction * Strain_energy_function_pt
Pointer to the strain energy function.
void calculate_second_piola_kirchhoff_stress(const DenseMatrix< double > &g, const DenseMatrix< double > &G, DenseMatrix< double > &sigma)
Calculate the contravariant 2nd Piola Kirchhoff stress tensor. Arguments are the covariant undeformed...
An OomphLibError object which should be thrown when an run-time error is encountered....
////////////////////////////////////////////////////////////////// //////////////////////////////////...
Definition: matrices.h:1701
virtual void derivatives(Vector< double > &I, Vector< double > &dWdI)
Return the derivatives of the strain energy function with respect to the strain invariants....
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const std::complex< double > I(0.0, 1.0)
The imaginary unit.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...