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-2025 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
31namespace 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 //========================================================================
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 //=======================================================================
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",
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(
93 }
94
95 // If the stress tensor does not have the same dimensions as the others
96 // complain.
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(
106 }
107 }
108
109
110 //===========================================================================
111 /// The function to calculate the contravariant tensor from a covariant one
112 //===========================================================================
115 {
116 // Initial error checking
117#ifdef PARANOID
118 // Test that the matrices are of the same dimension
120 {
121 throw OomphLibError("Matrices passed to calculate_contravariant() are "
122 "not of equal dimension",
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
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(
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",
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",
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 //===========================================================================
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
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(
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",
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",
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",
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,
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",
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
376
377 // Advanced metric tensor
378 DenseMatrix<double> G_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
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
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,
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",
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
479
480 // Advanced metric tensor etc.
481 DenseMatrix<double> G_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
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
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,
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",
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
593
594 // Advanced metric tensor etc
595 DenseMatrix<double> G_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
626
627 // Derivative of generalised dilatation
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
651 {
652 for (unsigned i = 0; i < dim; i++)
653 {
654 for (unsigned j = 0; j < i; j++)
655 {
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
697 // We don't need the Jacobian so cast the function to void
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,
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
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,
816 double& detG)
817 {
818 // Error checking
819#ifdef PARANOID
821#endif
822
823 // Find the dimension of the problem
824 unsigned dim = G.nrow();
825
826 // Calculate the contravariant Deformed metric tensor
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",
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);
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,
1014 double& detG)
1015 {
1016// Error checking
1017#ifdef PARANOID
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",
1035 }
1036#endif
1037
1038 // Calculate the contravariant undeformed and deformed metric tensors
1040 // Don't need this determinant
1042 // These are passed back
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,
1141 double& gen_dil,
1142 double& inv_kappa)
1143 {
1144// Error checking
1145#ifdef PARANOID
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",
1162 }
1163#endif
1164
1165 // Calculate the contravariant undeformed and deformed metric tensors
1166 // and get the determinants of the metric tensors
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...
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Definition elements.h:1185
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....
virtual void derivatives(Vector< double > &I, Vector< double > &dWdI)
Return the derivatives of the strain energy function with respect to the strain invariants....
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).