Merge pull request #6 from 3MAH/add_hyperelasticity

Add hyperelasticity constitutive models

include/simcoon/Continuum_mechanics/Functions/contimech.hpp

@@ -406,5 +406,40 @@ arma::mat auto_dyadic(const arma::mat &a);
 */
 arma::mat dyadic(const arma::mat &a, const arma::mat &b);
 
+/**
+ * @brief Provides the symmetric 4th-order dyadic product of a symmetric tensor with itself
+ * @param a
+ * @details This function returns the operation \f$ c = a \odot a \f$. The function returns a 6x6 matrix that correspond to a 4th order tensor. Note that such conversion to 6x6 matrices product correspond to a conversion with the component of the 4th order tensor correspond to the component of the matrix (such as stiffness matrices)
+\f[
+    \begin{align}
+    \left(\mathbf{a} \odot \mathbf{b} \right)_{ijkl} = \frac{1}{2} \left( a_{ik} a_{jl} + a_{il} a_{jk} \right)
+    \end{align}
+\f]
+ * @returns The 6x6 matrix that represent the dyadic product (arma::mat)
+ * @code 
+    mat a = randu(3,3);
+    mat b = randu(3,3);
+    mat c = sym_dyadic_operator(a,b);
+ * @endcode
+*/
+arma::mat auto_sym_dyadic_operator(const arma::mat &a);
+
+/**
+ * @brief Provides the symmetric 4th-order dyadic product of two symmetric tensors
+ * @param a, b
+ * @details This function returns the operation \f$ c = a \odot b \f$. The function returns a 6x6 matrix that correspond to a 4th order tensor. Note that such conversion to 6x6 matrices product correspond to a conversion with the component of the 4th order tensor correspond to the component of the matrix (such as stiffness matrices)
+\f[
+    \begin{align}
+    \left(\mathbf{a} \odot \mathbf{b} \right)_{ijkl} = \frac{1}{2} \left( a_{ik} b_{jl} + a_{il} b_{jk} \right)
+    \end{align}
+\f]
+ * @returns The 6x6 matrix that represent the dyadic product (arma::mat)
+ * @code 
+    mat a = randu(3,3);
+    mat b = randu(3,3);
+    mat c = sym_dyadic_operator(a,b);
+ * @endcode
+*/
+arma::mat sym_dyadic_operator(const arma::mat &a, const arma::mat &b);
 
 } //namespace simcoon

include/simcoon/Continuum_mechanics/Functions/objective_rates.hpp

@@ -501,7 +501,7 @@ arma::mat Dtau_LieDD_Dtau_JaumannDD(const arma::mat &Lt, const arma::mat &tau);
  * @brief Computes the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the logarithmic spin from the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis
  *
  * \f[
- *      \left( \frac{\partial \mathbf{\eth \tau}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} \mathbf{\Omega}_{\textrm{log}}} (i,s,p,r) = \left( \frac{\partial \mathbf{\eth \sigma}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\vec{g}_i)} (i,s,p,r) - \frac{1}{2} \tau (p,s) \delta (i,r) - \frac{1}{2} \tau (r,s) \delta (i,p) - \frac{1}{2} \tau (i,r) \delta (s,p) - \frac{1}{2} \tau(i,p) \delta(s,r)
+ *      \left( \frac{\partial \mathbf{\eth \tau}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} \mathbf{\Omega}_{\textrm{log}}} (i,s,p,r) = \left( \frac{\partial \mathbf{\eth \tau}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\vec{g}_i)} (i,s,p,r) - \frac{1}{2} \tau (p,s) \delta (i,r) - \frac{1}{2} \tau (r,s) \delta (i,p) - \frac{1}{2} \tau (i,r) \delta (s,p) - \frac{1}{2} \tau(i,p) \delta(s,r)
  * \f]
  * 
  * This function takes in the tangent modulus \f$ \left( \frac{\partial \mathbf{\eth \tau}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\vec{g}_i)} \f$ that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis,
@@ -518,24 +518,44 @@ arma::mat Dtau_LieDD_Dtau_JaumannDD(const arma::mat &Lt, const arma::mat &tau);
 arma::mat Dtau_LieDD_Dtau_logarithmicDD(const arma::mat &Lt, const arma::mat &BBBB, const arma::mat &tau);
 
 /**
- * @brief Computes the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin from the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis
+ * @brief Computes the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin from the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis
  *
  * \f[
- *      \left( \frac{\partial \mathbf{\eth \sigma}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\mathbf{W})} = \frac{1}{J} \left( \frac{\partial \mathbf{\eth \tau}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\mathbf{W})}
+ *      \left( \frac{\partial \mathbf{\eth \sigma}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\mathbf{W})} = \left( \frac{\partial \mathbf{\eth \sigma}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\mathbf{W})}
  * \f]
  * 
- * This function takes in the tangent modulus \f$ \left( \frac{\partial \mathbf{\eth \tau}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\vec{g}_i)} \f$ that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis,
+ * This function takes in the tangent modulus \f$ \left( \frac{\partial \mathbf{\eth \sigma}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\vec{g}_i)} \f$ that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis,
  * the logarithmic antisymmetric tensor-valued function \f$ \mathbf{\mathcal{B}}^{\textrm{log}} \f$, 
  * the transformation gradient \f$ \mathbf{F} \f$ and the Cauchy stress tensor \f$ \mathbf{\tau} \f$.
  * 
  * It returns the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin
  * 
  * @note : The return tensor is the tangent modulus utilized by Abaqus
  * 
- * @param[in] DSDE (6x6 arma::mat) tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$
- * @param[in] tau (3x3 arma::mat) Kirchoff stress tensor \f$ \mathbf{\tau} \f$.
+ * @param[in] Dsigma_LieDD (6x6 arma::mat) tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$
+ * @param[in] sigma (3x3 arma::mat) Kirchoff stress tensor \f$ \mathbf{\tau} \f$.
  * @return (6x6 arma::mat) the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the Zaremba-Jaumann-Noll spin
 */
-arma::mat Dtau_LieDD_Dsigma_JaumannDD(const arma::mat &Dtau_JaumannDD, const double &J);
-
+arma::mat Dsigma_LieDD_Dsigma_JaumannDD(const arma::mat &Dsigma_LieDD, const arma::mat &sigma);
+
+/**
+ * @brief Computes the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the logarithmic spin from the tangent modulus that links the Cauchy stress tensor \f$ \mathbf{\sigma} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis
+ *
+ * \f[
+ *      \left( \frac{\partial \mathbf{\eth \sigma}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} \mathbf{\Omega}_{\textrm{log}}} (i,s,p,r) = \left( \frac{\partial \mathbf{\eth \sigma}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\vec{g}_i)} (i,s,p,r) - \frac{1}{2} \sigma (p,s) \delta (i,r) - \frac{1}{2} \sigma (r,s) \delta (i,p) - \frac{1}{2} \tau (i,r) \delta (s,p) - \frac{1}{2} \tau(i,p) \delta(s,r)
+ * \f]
+ * 
+ * This function takes in the tangent modulus \f$ \left( \frac{\partial \mathbf{\eth \tau}}{\partial \eth \mathbf{e}} \right)_{\mathcal{R} (\vec{g}_i)} \f$ that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated in the natural covariant vector basis,
+ * the logarithmic antisymmetric tensor-valued function \f$ \mathbf{\mathcal{B}}^{\textrm{log}} \f$, 
+ * the transformation gradient \f$ \mathbf{F} \f$ and the Cauchy stress tensor \f$ \mathbf{\tau} \f$.
+ * 
+ * It returns the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the logarithmic spin
+ * 
+ * @param[in] DSDE (6x6 arma::mat) tangent modulus \f$ \frac{\partial \mathbf{S}}{\partial \mathbf{E}} \f$ that links the Piola-Kirchoff II stress \f$ \mathbf{S} \f$ to the Green-Lagrange stress \f$ \mathbf{E} \f$
+ * @param[in] BBBB (6x6 arma::mat) logarithmic antisymmetric tensor-valued function  \f$ \mathbf{\mathcal{B}}^{\textrm{log}} \f$
+ * @param[in] sigma (3x3 arma::mat) Kirchoff stress tensor \f$ \mathbf{\tau} \f$.
+ * @return (6x6 arma::mat) the tangent modulus that links the Kirchoff stress tensor \f$ \mathbf{\tau} \f$ and rate of deformation \f$ \mathbf{D} \f$ integrated using the logarithmic spin
+*/
+arma::mat Dsigma_LieDD_Dsigma_logarithmicDD(const arma::mat &Lt, const arma::mat &BBBB, const arma::mat &sigma);
+
 } //namespace simcoon

include/simcoon/Continuum_mechanics/Umat/Finite/generic_hyper_invariants.hpp

@@ -0,0 +1,37 @@
+/* This file is part of simcoon.
+ 
+ simcoon is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+ 
+ simcoon is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License
+ along with simcoon.  If not, see <http://www.gnu.org/licenses/>.
+ 
+ */
+
+///@file elastic_isotropic.hpp
+///@brief User subroutine for Isotropic elastic materials in 3D case
+///@version 1.0
+
+#pragma once
+
+#include <armadillo>
+
+namespace simcoon{
+
+///@brief The elastic UMAT requires 2 constants:
+///@brief props[0] : Young modulus
+///@brief props[1] : Poisson ratio
+///@brief props[2] : CTE
+
+///@brief No statev is required for thermoelastic constitutive law
+
+void umat_generic_hyper_invariants(const std::string &umat_name, const arma::vec &, const arma::vec &, const arma::mat &, const arma::mat&, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &);
+
+} //namespace simcoon

include/simcoon/Continuum_mechanics/Umat/Finite/generic_hyper.hpp → include/simcoon/Continuum_mechanics/Umat/Finite/generic_hyper_pstretch.hpp

@@ -32,6 +32,6 @@ namespace simcoon{
 
 ///@brief No statev is required for thermoelastic constitutive law
 
-void umat_generic_hyper(const arma::vec &, const arma::vec &, const arma::mat &, const arma::mat&, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &);
+void umat_generic_hyper_pstretch(const arma::vec &, const arma::vec &, const arma::mat &, const arma::mat&, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &);
 
 } //namespace simcoon

src/Continuum_mechanics/Functions/contimech.cpp

@@ -409,4 +409,79 @@ mat dyadic(const mat &A, const mat &B) {
     return FTensor4_mat(C_);    
 }
 
+///This computes the symmetric 4th-order dyadic product A o A = 0.5*(A(i,j)*A(k,l) + A(i,j)*A(l,k));
+mat auto_sym_dyadic_operator(const mat &A) {
+
+	mat C = zeros(6,6);
+
+	int ij=0;
+	int kl=0;
+
+    umat Id(3,3);
+    Id(0,0) = 0;
+    Id(0,1) = 3;
+    Id(0,2) = 4;
+    Id(1,0) = 3;
+    Id(1,1) = 1;
+    Id(1,2) = 5;
+    Id(2,0) = 4;
+    Id(2,1) = 5;
+    Id(2,2) = 2;
+
+	for (int i=0; i<3; i++) {
+		for (int j=i; j<3; j++) {
+			ij = Id(i,j);
+			for (int k=0; k<3; k++) {
+				for (int l=k; l<3; l++) {
+					kl = Id(k,l);
+					C(ij,kl) = 0.5*(A(i,k)*A(j,l) + A(i,l)*A(j,k));
+				}
+			}
+		}
+	}
+
+	return C;
+}
+
+///This computes the symmetric 4th-order dyadic product A o B = 0.5*(A(i,j)*B(k,l) + A(i,j)*B(l,k));
+mat sym_dyadic_operator(const mat &A, const mat &B) {
+
+	mat C = zeros(6,6);
+
+	int ij=0;
+	int kl=0;
+
+    umat Id(3,3);
+    Id(0,0) = 0;
+    Id(0,1) = 3;
+    Id(0,2) = 4;
+    Id(1,0) = 3;
+    Id(1,1) = 1;
+    Id(1,2) = 5;
+    Id(2,0) = 4;
+    Id(2,1) = 5;
+    Id(2,2) = 2;
+
+	for (int i=0; i<3; i++) {
+		for (int j=i; j<3; j++) {
+			ij = Id(i,j);
+			for (int k=0; k<3; k++) {
+				for (int l=k; l<3; l++) {
+					kl = Id(k,l);
+					C(ij,kl) = 0.5*(A(i,k)*B(j,l) + A(i,l)*B(j,k));
+				}
+			}
+		}
+	}
+	return C;
+}
+
+/*mat eulerian_determinant(const mat &A) {
+	mat Id = eye(3,3);
+	vec A_Id = A*Ith();
+	double Id_A_Id = sum(Ith()%A_Id);
+	return A - (1./3.)*(sym_dyadic(Id,A_Id)+sym_dyadic(A_Id,Id))+(1./9.)*Id_A_Id*auto_sym_dyadic(Id);
+
+}*/
+
 } //namespace simcoon