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