/*
 * Fibre Dispersion Law
 * 
 * Model Structure
 * 
 * Constitutive relations are fundamental to the solution of problems
 * in continuum mechanics, and are required in the study of, for
 * example, mechanically dominated clinical interventions involving
 * soft biological tissues. Structural continuum constitutive models
 * of arterial layers integrate information about the tissue morphology
 * and therefore allow investigation of the interrelation between
 * structure and function in response to mechanical loading. Collagen
 * fibres are key ingredients in the structure of arteries. In
 * the media (the middle layer of the artery wall) they are arranged
 * in two helically distributed families with a small pitch and
 * very little dispersion in their orientation (i.e. they are aligned
 * quite close to the circumferential direction). By contrast,
 * in the adventitial and intimal layers, the orientation of the
 * collagen fibres is dispersed, as shown by polarized light microscopy
 * of stained arterial tissue. As a result, continuum models that
 * do not account for the dispersion are not able to capture accurately
 * the stress-strain response of these layers. The purpose of this
 * paper, therefore, is to develop a structural continuum framework
 * that is able to represent the dispersion of the collagen fibre
 * orientation. This then allows the development of a new hyperelastic
 * free-energy function that is particularly suited for representing
 * the anisotropic elastic properties of adventitial and intimal
 * layers of arterial walls, and is a generalization of the fibre-reinforced
 * structural model introduced by Holzapfel and Gasser (Holzapfel
 * and Gasser 2001 Comput. Meth. Appl. Mech. Eng. 190, 4379-4403)
 * and Holzapfel et al. (Holzapfel et al. 2000 J. Elast. 61, 1-48).
 * The model incorporates an additional scalar structure parameter
 * that characterizes the dispersed collagen orientation. An efficient
 * finite element implementation of the model is then presented
 * and numerical examples show that the dispersion of the orientation
 * of collagen fibres in the adventitia of human iliac arteries
 * has a significant effect on their mechanical response.
 * 
 * The original paper reference is cited below:
 * 
 * Hyperelastic modelling of arterial layers with distributed collagen
 * fibre orientations, T.C. Gasser, R.W. Ogden and G.A. Holzapfel,
 * 2006, Journal of the Royal Society, 3(6), 15-35. PubMed ID:
 * 16849214
 */

import nsrunit;
unit conversion on;
unit strain=1  dimensionless;
unit stress=1  dimensionless;

math main {
	//Warning:  the following variables were set 'extern' or given
	//  an initial value of '0' because the model would otherwise be
	//  underdetermined:  Emm, Emn, Ems, Enn, Ens, Ess
	extern real Emm strain;
	extern real Emn strain;
	extern real Ems strain;
	extern real Enn strain;
	extern real Ens strain;
	extern real Ess strain;
	real c11 strain;
	c11=0;
	real c12 strain;
	c12=0;
	real c21 strain;
	c21=0;
	real c22 strain;
	c22=0;
	real c41 strain;
	c41=0;
	real c42 strain;
	c42=0;
	real c81 strain;
	c81=0;
	real c82 strain;
	c82=0;
	real kappa1 strain;
	kappa1=0;
	real kappa2 strain;
	kappa2=0;
	real Tdevmm stress;
	real Tdevmn stress;
	real Tdevms stress;
	real Tdevnn stress;
	real Tdevns stress;
	real Tdevss stress;
	real I1 strain;
	real I2 strain;
	real I4 strain;
	real I6 strain;
	real I8 strain;

	// <component name="interface">

	// <component name="equations">
	I1=(3+2*Emm+2*Enn+2*Ess);
	I2=(3+(-4)*Emn*Emn+4*Enn+(-4)*Ems*Ems+(-4)*Ens*Ens+4*(1+Enn)*Ess+4*(1+Enn)*Ess+4*(1+Enn+Ess)*Emm);
	I4=(1+2*Emm);
	I6=(1+2*Enn);
	I8=(1+2*Ess);
	Tdevmm=(2*c11*exp(c12*(I1-3)^2)*(I1-3)+4*c21*exp(c22*(I2-3)^2)*(I2-3)*(1+Enn+Ess)+2*c41*exp(c42*(I4*(1-3*kappa1)-1)^2)*(I4*(1-3*kappa1)-1)*(1-3*kappa1));
	Tdevmn=((-4)*c21*exp(c22*(I2-3)^2)*(I2-3)*Emn);
	Tdevms=((-4)*c21*exp(c22*(I2-3)^2)*(I2-3)*Ems);
	Tdevnn=(2*c11*exp(c12*(I1-3)^2)*(I1-3)+4*c21*exp(c22*(I2-3)^2)*(I2-3)*(1+Emm+Ess));
	Tdevss=(2*c11*exp(c12*(I1-3)^2)*(I1-3)+4*c21*exp(c22*(I2-3)^2)*(I2-3)*(1+Emm+Enn)+2*c81*exp(c82*(I8*(1-3*kappa2)-1)^2)*(I8*(1-3*kappa2)-1)*(1-3*kappa2));
	Tdevns=((-4)*c21*exp(c22*(I2-3)^2)*(I2-3)*Ens);
}