39 if (M.nrow() != M.ncol())
57 if ((
M1.nrow() !=
M2.nrow()) || (
M1.ncol() !=
M2.ncol()))
87 std::string error_message =
"Deformed metric tensor does \n";
89 "not have same dimensions as the undeformed metric tensor\n";
99 std::string error_message =
100 "Strain tensor passed to calculate_green_strain() does \n";
102 "not have same dimensions as the undeformed metric tensor\n";
121 throw OomphLibError(
"Matrices passed to calculate_contravariant() are "
122 "not of equal dimension",
129 unsigned dim =
Gdown.ncol();
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";
209 throw OomphLibError(
"Dimension of matrix must be 0, 1, 2 or 3\n",
231 const unsigned dim =
Gdown.ncol();
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";
309 "Analytic derivatives of metric tensors not yet implemented in 3D\n",
334 throw OomphLibError(
"Dimension of matrix must be 0, 1, 2 or 3\n",
364 throw OomphLibError(
"Matrices passed are not of equal dimension",
371 const unsigned dim =
G.ncol();
382 for (
unsigned i = 0;
i < dim;
i++)
384 for (
unsigned j = 0;
j < dim;
j++)
396 for (
unsigned i = 0;
i < dim;
i++)
398 for (
unsigned j =
i;
j < dim;
j++)
406 for (
unsigned ii = 0;
ii < dim;
ii++)
408 for (
unsigned jj =
ii;
jj < dim;
jj++)
424 for (
unsigned i = 0;
i < dim;
i++)
426 for (
unsigned j = 0;
j <
i;
j++)
428 for (
unsigned ii = 0;
ii < dim;
ii++)
430 for (
unsigned jj = 0;
jj <
ii;
jj++)
458 const double& interpolated_solid_p,
468 throw OomphLibError(
"Matrices passed are not of equal dimension",
475 const unsigned dim =
G.ncol();
487 for (
unsigned i = 0;
i < dim;
i++)
489 for (
unsigned j = 0;
j < dim;
j++)
501 for (
unsigned i = 0;
i < dim;
i++)
503 for (
unsigned j =
i;
j < dim;
j++)
518 for (
unsigned ii = 0;
ii < dim;
ii++)
520 for (
unsigned jj =
ii;
jj < dim;
jj++)
538 for (
unsigned i = 0;
i < dim;
i++)
540 for (
unsigned j = 0;
j <
i;
j++)
544 for (
unsigned ii = 0;
ii < dim;
ii++)
546 for (
unsigned jj = 0;
jj <
ii;
jj++)
572 const double& interpolated_solid_p,
582 throw OomphLibError(
"Matrices passed are not of equal dimension",
589 const unsigned dim =
G.ncol();
602 for (
unsigned i = 0;
i < dim;
i++)
604 for (
unsigned j = 0;
j < dim;
j++)
616 for (
unsigned i = 0;
i < dim;
i++)
618 for (
unsigned j =
i;
j < dim;
j++)
632 for (
unsigned ii = 0;
ii < dim;
ii++)
634 for (
unsigned jj =
ii;
jj < dim;
jj++)
652 for (
unsigned i = 0;
i < dim;
i++)
654 for (
unsigned j = 0;
j <
i;
j++)
658 for (
unsigned ii = 0;
ii < dim;
ii++)
660 for (
unsigned jj = 0;
jj <
ii;
jj++)
693 unsigned dim =
G.nrow();
701 double C1 = (*E_pt) / (2.0 * (1.0 + (*Nu_pt)));
702 double C2 = 2.0 * (*Nu_pt) / (1.0 - 2.0 * (*Nu_pt));
708 for (
unsigned i = 0;
i < dim;
i++)
710 for (
unsigned j =
i;
j < dim;
j++)
717 for (
unsigned i = 0;
i < dim;
i++)
719 for (
unsigned j = 0;
j <
i;
j++)
726 for (
unsigned i = 0;
i < dim;
i++)
728 for (
unsigned j =
i;
j < dim;
j++)
732 for (
unsigned k = 0;
k < dim;
k++)
734 for (
unsigned l = 0;
l < dim;
l++)
746 for (
unsigned i = 0;
i < dim;
i++)
748 for (
unsigned j = 0;
j <
i;
j++)
750 sigma(
i,
j) = sigma(
j,
i);
774 unsigned dim =
G.nrow();
785 (1.0 - 2.0 * (*Nu_pt)) * (1.0 + (*
Nu_pt)) / ((*
E_pt) * (*Nu_pt));
790 for (
unsigned i = 0;
i < dim;
i++)
792 for (
unsigned j = 0;
j < dim;
j++)
824 unsigned dim =
G.nrow();
830 double C1 = (*E_pt) / (2.0 * (1.0 + (*Nu_pt)));
836 for (
unsigned i = 0;
i < dim;
i++)
838 for (
unsigned j =
i;
j < dim;
j++)
845 for (
unsigned i = 0;
i < dim;
i++)
847 for (
unsigned j = 0;
j <
i;
j++)
854 for (
unsigned i = 0;
i < dim;
i++)
856 for (
unsigned j =
i;
j < dim;
j++)
860 for (
unsigned k = 0;
k < dim;
k++)
862 for (
unsigned l = 0;
l < dim;
l++)
873 for (
unsigned i = 0;
i < dim;
i++)
875 for (
unsigned j = 0;
j <
i;
j++)
906 unsigned dim = g.nrow();
912 "IsotropicStrainEnergyFunctionConstitutiveLaw::";
913 function_name +=
"calculate_second_piola_kirchhoff_stress()";
915 throw OomphLibError(
"Check constitutive equations carefully when dim=1",
932 for (
unsigned i = 0;
i < dim;
i++)
934 for (
unsigned j = 0;
j < dim;
j++)
964 if (std::fabs(
dWdI[1]) > 0.0)
966 for (
unsigned i = 0;
i < dim;
i++)
968 for (
unsigned j = 0;
j < dim;
j++)
971 for (
unsigned r = 0;
r < dim;
r++)
973 for (
unsigned s = 0;
s < dim;
s++)
988 double p = 2.0 *
dWdI[2] * I[2];
991 for (
unsigned i = 0;
i < dim;
i++)
993 for (
unsigned j = 0;
j < dim;
j++)
1022 unsigned dim = g.nrow();
1029 "IsotropicStrainEnergyFunctionConstitutiveLaw::";
1030 function_name +=
"calculate_second_piola_kirchhoff_stress()";
1032 throw OomphLibError(
"Check constitutive equations carefully when dim=1",
1050 for (
unsigned i = 0;
i < dim;
i++)
1052 for (
unsigned j = 0;
j < dim;
j++)
1078 if (std::fabs(
dWdI[1]) > 0.0)
1080 for (
unsigned i = 0;
i < dim;
i++)
1082 for (
unsigned j = 0;
j < dim;
j++)
1085 for (
unsigned r = 0;
r < dim;
r++)
1087 for (
unsigned s = 0;
s < dim;
s++)
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)));
1113 K = (I[0] *
phi + 2.0 * I[1] *
psi) / 3.0;
1119 for (
unsigned i = 0;
i < dim;
i++)
1121 for (
unsigned j = 0;
j < dim;
j++)
1150 unsigned dim = g.nrow();
1156 "IsotropicStrainEnergyFunctionConstitutiveLaw::";
1157 function_name +=
"calculate_second_piola_kirchhoff_stress()";
1159 throw OomphLibError(
"Check constitutive equations carefully when dim=1",
1176 for (
unsigned i = 0;
i < dim;
i++)
1178 for (
unsigned j = 0;
j < dim;
j++)
1207 if (std::fabs(
dWdI[1]) > 0.0)
1209 for (
unsigned i = 0;
i < dim;
i++)
1211 for (
unsigned j = 0;
j < dim;
j++)
1214 for (
unsigned r = 0;
r < dim;
r++)
1216 for (
unsigned s = 0;
s < dim;
s++)
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)));
1244 K = (I[0] *
phi + 2.0 * I[1] *
psi) / 3.0;
1258 for (
unsigned i = 0;
i < dim;
i++)
1260 for (
unsigned j = 0;
j < dim;
j++)
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.
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).