{ "cells": [ { "cell_type": "markdown", "id": "5069db8e-28d6-4aec-a8c8-1a46b224548d", "metadata": {}, "source": [ "# Using MatModel" ] }, { "cell_type": "code", "execution_count": 1, "id": "9f45ba17-b251-4ec8-ac1a-2e3f493076e5", "metadata": {}, "outputs": [], "source": [ "import pymecht as pmt\n", "import numpy as np" ] }, { "cell_type": "markdown", "id": "2b869706-e84f-424a-87a2-ad05dc5595e6", "metadata": {}, "source": [ "## Single material model" ] }, { "cell_type": "markdown", "id": "60fd3f30-4489-4a75-b483-d38852f47413", "metadata": {}, "source": [ "The core module in pyMechT is `MatModel` which allows us to add models together to simulate the stress-strain behavior. A basic usage is when there is only one model, for example we can create a material with Yeoh model, which is isotropic and only a function of the first invariant:" ] }, { "cell_type": "code", "execution_count": 2, "id": "b1f8657e-7ede-4c82-8b33-d0430bf8dc35", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "------------------------------------------------------------------\n", "Keys Value Fixed? Lower bound Upper bound \n", "------------------------------------------------------------------\n", "c1_0 1.00 No 1.00e-04 1.00e+02 \n", "c2_0 1.00 No 0.00 1.00e+02 \n", "c3_0 1.00 No 0.00 1.00e+02 \n", "c4_0 0.00 No 0.00 1.00e+02 \n", "------------------------------------------------------------------\n", "\n" ] } ], "source": [ "mat1 = pmt.MatModel('yeoh')\n", "print(mat1.parameters) # Returns parameters as a dictionary" ] }, { "cell_type": "markdown", "id": "83025e51-0835-49f8-b415-7ee060729bae", "metadata": {}, "source": [ "In principle, one can directly use this material to calculate stresses given a deformation gradient using `stress` method" ] }, { "cell_type": "code", "execution_count": 3, "id": "bfecaddb-a27d-4ed8-9b7a-e6cafa57cd29", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mSignature:\u001b[0m\n", "\u001b[0mmat1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstress\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mF\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mtheta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mstresstype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'cauchy'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mincomp\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m \u001b[0mFdiag\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", "\u001b[0;34m\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mDocstring:\u001b[0m\n", "Returns the stress tensor of the material model\n", "\n", "Parameters\n", "----------\n", "F: np.array([3,3])\n", " the deformation gradient, default: identity matrix (no deformation)\n", "\n", "theta: ParamDict or dict\n", " the parameters of the model, if None, then the default values are used\n", "\n", "stresstype: str\n", " the type of stress tensor to return with the following options (case insensitive) \n", "\n", " * 'cauchy': Cauchy stress,\n", " * '1pk' or '1stpk' or 'firstpk': 1st Piola-Kirchoff stress\n", " * '2pk' or '2ndpk' or 'secondpk': 2nd Piola-Kirchoff stress\n", " * default: Cauchy stress\n", " \n", "incomp: bool\n", " if True, then the material is assumed to be incompressible, default: False\n", " \n", "Fdiag: bool\n", " if True, then it is assumed that F is diagonal (for faster computation), default: False\n", "\n", "Returns\n", "-------\n", "np.array([3,3])\n", " the stress 3X3 tensor\n", "\u001b[0;31mFile:\u001b[0m /usr/local/lib/python3.9/site-packages/pymecht/MatModel.py\n", "\u001b[0;31mType:\u001b[0m method\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mat1.stress?" ] }, { "cell_type": "markdown", "id": "3db154a7-d75b-43f8-8d0c-e345cd4a0733", "metadata": {}, "source": [ "Thus, we can calculate the Cauchy stress at an identity deformation gradient. However, since this model does not have a volumetric part, we must set incompressibility to be true to get a zero stress. This is done using a Lagrange multiplier, which is calculated internally in pyMechT by setting the normal stress along the third axis equal to zero." ] }, { "cell_type": "code", "execution_count": 4, "id": "b0e3992d-f494-4cdd-8ca2-f631dd58b658", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[0., 0., 0.],\n", " [0., 0., 0.],\n", " [0., 0., 0.]])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mat1.stress(np.eye(3),incomp=True)" ] }, { "cell_type": "markdown", "id": "b5bfc716-1340-44c9-9ffa-ed58fa4eeed7", "metadata": {}, "source": [ "Or, we can calculate other stresses at non-identity deformation gradient" ] }, { "cell_type": "code", "execution_count": 5, "id": "e2eb28a0-d013-4768-aebd-8b7e5b3951be", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[3.78128667 0.71615278 0. ]\n", " [1.03126 3.78128667 0. ]\n", " [0. 0. 0. ]]\n", "[[3.15107222 0.59679398 0. ]\n", " [0.59679398 3.10133939 0. ]\n", " [0. 0. 0. ]]\n" ] } ], "source": [ "F = np.array([[1.2,0,0],[0.1,1.2,0],[0,0,1]])\n", "print(mat1.stress(F,stresstype='1pk',incomp=True))\n", "print(mat1.stress(F,stresstype='2pk',incomp=True))" ] }, { "cell_type": "markdown", "id": "d4129361-dd57-46df-8d9c-dc3cecacdb95", "metadata": {}, "source": [ "One can also change the parameters and pass them while calculating the stress." ] }, { "cell_type": "code", "execution_count": 6, "id": "c563a34d-4a39-45ad-90a4-1dca1350bbc4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[4.51462 0.85504167 0. ]\n", " [1.23126 4.51462 0. ]\n", " [0. 0. 0. ]]\n" ] } ], "source": [ "params = mat1.parameters\n", "params['c1_0'].set(2)\n", "print(mat1.stress(F,params,stresstype='1pk',incomp=True))" ] }, { "cell_type": "markdown", "id": "e6858cf7-fb2b-4b13-8e33-9c12e0971cce", "metadata": {}, "source": [ "Alternatively, instead of passing the parameters, one can set the values and call the `stress` function without the parameters" ] }, { "cell_type": "code", "execution_count": 7, "id": "2397eb8c-010e-4352-adb8-8e2bf38030bb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[4.51462 0.85504167 0. ]\n", " [1.23126 4.51462 0. ]\n", " [0. 0. 0. ]]\n" ] } ], "source": [ "mat1.parameters = params\n", "print(mat1.stress(F,stresstype='1pk',incomp=True))" ] }, { "cell_type": "markdown", "id": "82274a7b-6bdf-4c3b-9788-2d1c64bb9e0b", "metadata": {}, "source": [ "## Adding materials" ] }, { "cell_type": "markdown", "id": "d698c929-c1b5-4bbe-b7c1-316cbb9114de", "metadata": {}, "source": [ "A convenient feature of pyMechT is that we can easily add different models together (potentially, the same model repeated). For example" ] }, { "cell_type": "code", "execution_count": 8, "id": "5173b326-fd60-446b-b582-2f682d6b984c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Material model with 2 components:\n", "Component1: YEOH\n", "Component2: NH\n", " ------------------------------------------------------------------\n", "Keys Value Fixed? Lower bound Upper bound \n", "------------------------------------------------------------------\n", "c1_0 1.00 No 1.00e-04 1.00e+02 \n", "c2_0 1.00 No 0.00 1.00e+02 \n", "c3_0 1.00 No 0.00 1.00e+02 \n", "c4_0 0.00 No 0.00 1.00e+02 \n", "mu_1 1.00 No 1.00e-04 1.00e+02 \n", "------------------------------------------------------------------\n", "\n", "Material model with 3 components:\n", "Component1: GOH\n", "Component2: GOH\n", "Component3: NH\n", " ------------------------------------------------------------------\n", "Keys Value Fixed? Lower bound Upper bound \n", "------------------------------------------------------------------\n", "k1_0 10.00 No 0.10 30.00 \n", "k2_0 10.00 No 0.10 30.00 \n", "k3_0 0.10 No 0.00 0.33 \n", "k1_1 10.00 No 0.10 30.00 \n", "k2_1 10.00 No 0.10 30.00 \n", "k3_1 0.10 No 0.00 0.33 \n", "mu_2 1.00 No 1.00e-04 1.00e+02 \n", "------------------------------------------------------------------\n", "\n" ] } ], "source": [ "mat1 = pmt.MatModel('yeoh','nh')\n", "mat2 = pmt.MatModel('goh','goh','nh')\n", "print(mat1, mat1.parameters)\n", "print(mat2, mat2.parameters)" ] }, { "cell_type": "markdown", "id": "db102528-0d67-48b6-bc7e-dbe107d17ccb", "metadata": {}, "source": [ "We can even add them afterwards, although this is rarely needed" ] }, { "cell_type": "code", "execution_count": 9, "id": "47d67913-4acd-43bd-96de-6425fd1e93d6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Material model with 5 components:\n", "Component1: YEOH\n", "Component2: NH\n", "Component3: GOH\n", "Component4: GOH\n", "Component5: NH\n", " ------------------------------------------------------------------\n", "Keys Value Fixed? Lower bound Upper bound \n", "------------------------------------------------------------------\n", "c1_0 1.00 No 1.00e-04 1.00e+02 \n", "c2_0 1.00 No 0.00 1.00e+02 \n", "c3_0 1.00 No 0.00 1.00e+02 \n", "c4_0 0.00 No 0.00 1.00e+02 \n", "mu_1 1.00 No 1.00e-04 1.00e+02 \n", "k1_2 10.00 No 0.10 30.00 \n", "k2_2 10.00 No 0.10 30.00 \n", "k3_2 0.10 No 0.00 0.33 \n", "k1_3 10.00 No 0.10 30.00 \n", "k2_3 10.00 No 0.10 30.00 \n", "k3_3 0.10 No 0.00 0.33 \n", "mu_4 1.00 No 1.00e-04 1.00e+02 \n", "------------------------------------------------------------------\n", "\n" ] } ], "source": [ "mat3 = mat1 + mat2\n", "print(mat3,mat3.parameters)" ] }, { "cell_type": "markdown", "id": "5d65d317-201a-4ead-aefd-abd1b5dee492", "metadata": {}, "source": [ "To get the model components back, we can use the `.models` " ] }, { "cell_type": "code", "execution_count": 10, "id": "0f71bf25-fa6a-4ee8-8a0e-8f639d05e6d9", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(,\n", " ,\n", " )" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mat2.models" ] }, { "cell_type": "markdown", "id": "f039fb65-47cd-46b5-9486-b0bd470a37d0", "metadata": {}, "source": [ "## Setting fiber directions" ] }, { "cell_type": "markdown", "id": "55e73d80-b01c-40b6-aa8a-aa2585373618", "metadata": {}, "source": [ "Note that, for anisotropic models (such as `goh`), fiber directions need to be specified before they can be called. That is, `mat1` above can be used but not `mat2`" ] }, { "cell_type": "code", "execution_count": 11, "id": "cc52dd71-3285-4cce-854e-339a33d865fb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0. 0. 0.]\n", " [0. 0. 0.]\n", " [0. 0. 0.]]\n", "Value error occured GOH model class uses I4 but no fiber directions have been defined. Did you forget to set the fiber directions?\n" ] } ], "source": [ "print(mat1.stress(np.eye(3),incomp=True))\n", "try:\n", " mat2.stress(np.eye(3),incomp=True)\n", "except ValueError as e:\n", " print(\"Value error occured\", e)" ] }, { "cell_type": "markdown", "id": "cf9e40be-de45-4810-9a5d-4b50c6cd7ab8", "metadata": {}, "source": [ "The fiber directions are set for each component model, and each component can have multiple fiber directions (the response is summed over them). This means, that the individual components can have different fiber directions (which is why one might want to add two GOH models together). Below we set two fiber directions for each of the GOH component." ] }, { "cell_type": "code", "execution_count": 12, "id": "d6316207-a53e-4b96-8d83-926c591124bd", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Material model with 3 components:\n", "Component1: GOH with fiber direction(s):[array([0.70710678, 0.70710678, 0. ]), array([0., 1., 0.])]\n", "Component2: GOH with fiber direction(s):[array([1., 0., 0.]), array([0.4472136 , 0.89442719, 0. ])]\n", "Component3: NH\n", "\n", "[[0. 0. 0.]\n", " [0. 0. 0.]\n", " [0. 0. 0.]]\n", "[[ 84.11760455 66.84260966 0. ]\n", " [ 66.84260966 126.67097888 0. ]\n", " [ 0. 0. 0. ]]\n" ] } ], "source": [ "model_comps = mat2.models\n", "model_comps[0].fiber_dirs = [ np.array([1,1,0]), np.array([0,1,0])]\n", "model_comps[1].fiber_dirs = [ np.array([1,0,0]), np.array([0.5,1,0])]\n", "print(mat2)\n", "print(mat2.stress(np.eye(3),incomp=True))\n", "print(mat2.stress(F,incomp=True))" ] }, { "cell_type": "markdown", "id": "a024c019-9df3-428f-89fc-16e5e2bbcba9", "metadata": {}, "source": [ "Note that the fiber direction vectors we supplied were not unit vectors. Internally they are made unit vectors. Thus, when we print them, we see some differences. Lastly, if we also set fiber directions to isotropic materials, it does not affect their response." ] }, { "cell_type": "code", "execution_count": 13, "id": "677da234-a423-410b-b0b7-3d84fec50b9d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Material model with 3 components:\n", "Component1: GOH with fiber direction(s):[array([0.70710678, 0.70710678, 0. ]), array([0., 1., 0.])]\n", "Component2: GOH with fiber direction(s):[array([1., 0., 0.]), array([0.4472136 , 0.89442719, 0. ])]\n", "Component3: NH with fiber direction(s):[array([1., 0., 0.]), array([0.4472136 , 0.89442719, 0. ])]\n", "\n", "[[0. 0. 0.]\n", " [0. 0. 0.]\n", " [0. 0. 0.]]\n", "[[ 84.11760455 66.84260966 0. ]\n", " [ 66.84260966 126.67097888 0. ]\n", " [ 0. 0. 0. ]]\n" ] } ], "source": [ "model_comps[2].fiber_dirs = [ np.array([1,0,0]), np.array([0.5,1,0])]\n", "print(mat2)\n", "print(mat2.stress(np.eye(3),incomp=True))\n", "print(mat2.stress(F,incomp=True))" ] }, { "cell_type": "markdown", "id": "fef2f57a-4a9d-414d-971f-31a9c143192c", "metadata": {}, "source": [ "However, for using pyMechT, one would seldom use the `MatModel` directly. Instead, it would be simply created and then used to create a `SampleExperiment` instance, which internally sets the deformation gradient based on the mode of deformation, incompressibility condition, and even a helper function to set the fiber directions more easily. These aspects are covered in the next examples." ] }, { "cell_type": "markdown", "id": "c692afab-b06b-4f02-a059-e20a1cf449eb", "metadata": {}, "source": [ "## Arbitrary `ARB` model" ] }, { "cell_type": "markdown", "id": "3b766da5-799c-4e16-8feb-d099028b4548", "metadata": {}, "source": [ "If the desired hyperelastic model is not implemented in `pymecht`, then one can use the arbitrary `ARB` model. It uses symbolic toolbox `sympy` to convert a string in terms of `I1`, `I2`, `J` and `I4` into a material model. An example is follows." ] }, { "cell_type": "code", "execution_count": 14, "id": "00137c16-4618-4ab2-9254-319212aaf372", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mInit signature:\u001b[0m \u001b[0mpmt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mARB\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_W\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_init_guess\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_low_bound\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_up_bound\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdparams\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mDocstring:\u001b[0m \n", "A material model class that allows the arbitrary definition of any strain energy density function (SEDF).\n", "Sympy's symbolic differentiation is used to calculate the partial derivatives of the SEDF with respect to the invariants I1, I2, J, and I4.\n", "The SEDF is provided as a string as a function of the invariants I1, I2, J, and I4, alongside initial guesses and upper/lower bounds for the parameters.\n", "The SEDF an initial guess for the parameters must be provided to identify parameters of the material model.\n", "If the strings are not provided, the user will be prompted to provide them.\n", "\n", "The parsed strings have the following functions replaced by their respective numpy equivalents:\n", " * exp -> np.exp\n", " * sqrt -> np.sqrt\n", " * log -> np.log\n", " * log10 -> np.log10\n", " * log2 -> np.log2\n", " * sin -> np.sin\n", " * cos -> np.cos\n", " * tan -> np.tan\n", " * arcsin -> np.arcsin\n", " * asin -> np.asin\n", " * arccos -> np.arccos\n", " * acos -> np.acos\n", " * arctan -> np.arctan\n", " * atan -> np.atan\n", " * hypot -> np.hypot\n", " * arctan2 -> np.arctan2\n", " * sinh -> np.sinh\n", " * cosh -> np.cosh\n", " * tanh -> np.tanh\n", " * arcsinh -> np.arcsinh\n", " * asinh -> np.asinh\n", " * arccosh -> np.arccosh\n", " * acosh -> np.acosh\n", " * arctanh -> np.arctanh\n", " * atanh -> np.atanh\n", "\n", "Example\n", "-------\n", " >>> from MatModel import *\n", " >>> model = pmt.ARB('mu/2.*(I1-3)','mu=1.','mu=0.01','mu=10.')\n", " >>> mat = pmt.MatModel(model)\n", " >>> F = np.random.rand\n", " >>> mat.stress(F)\n", "\u001b[0;31mFile:\u001b[0m /usr/local/lib/python3.9/site-packages/pymecht/MatModel.py\n", "\u001b[0;31mType:\u001b[0m type\n", "\u001b[0;31mSubclasses:\u001b[0m \n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pmt.ARB?" ] }, { "cell_type": "code", "execution_count": 15, "id": "638e3c73-c02d-4b95-a189-c936b1b7ed8e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Material model with 1 component:\n", "Component1: ARB\n", " ------------------------------------------------------------------\n", "Keys Value Fixed? Lower bound Upper bound \n", "------------------------------------------------------------------\n", "mu_0 1.00 No 1.00e-02 10.00 \n", "nu_0 1.00 No 1.00e-02 10.00 \n", "------------------------------------------------------------------\n", "\n", "[[0. 0. 0.]\n", " [0. 0. 0.]\n", " [0. 0. 0.]]\n", "[[2.27637222 0.55893333 0. ]\n", " [0.55893333 2.32295 0. ]\n", " [0. 0. 0. ]]\n" ] } ], "source": [ "model_I1I2 = pmt.ARB('mu/2.*(I1-3.)+nu*(I1-3)*(I2-3)','mu=1.,nu=1.','mu=0.01, nu=0.01','mu=10., nu=10.')\n", "mat = pmt.MatModel(model_I1I2)\n", "print(mat, mat.parameters)\n", "print(mat.stress(np.eye(3),incomp=True))\n", "print(mat.stress(F,incomp=True))" ] }, { "cell_type": "markdown", "id": "8d1f770e-a3de-4595-a678-f8014b3764e3", "metadata": {}, "source": [ "## Spline based data-driven model" ] }, { "cell_type": "markdown", "id": "cf520b7a-6f5f-482d-ab98-9a252ce730d9", "metadata": {}, "source": [ "Instead of an analytical expression, if one has data of the strain energy density function, it can be interpolated with a spline and used as a material model. `pymecht` has `splineI1` (a function of I1 alone) and `splineI1I4` (a function of I1 and I4) implemented to achieve this. These can be used as follows:" ] }, { "cell_type": "code", "execution_count": 16, "id": "9b7f34ac-c4e5-40d4-8a65-0b6c17a089bb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Material model with 1 component:\n", "Component1: splineI1\n", " ------------------------------------------------------------------\n", "Keys Value Fixed? Lower bound Upper bound \n", "------------------------------------------------------------------\n", "alpha_0 1.00 No -1.00e+01 10.00 \n", "------------------------------------------------------------------\n", "\n", "[[1.10778312 0.30212267 0. ]\n", " [0.30212267 1.13296001 0. ]\n", " [0. 0. 0. ]]\n" ] } ], "source": [ "from scipy.interpolate import make_interp_spline, RectBivariateSpline\n", "\n", "I1 = np.linspace(3,3.1,10)\n", "Psi = np.tan(I1-3)\n", "sp = make_interp_spline(I1,Psi) #create an interpolating spline\n", "mat = pmt.MatModel('splineI1')\n", "mat.models[0].set(sp) #set the spline for the model\n", "print(mat,mat.parameters)\n", "print(mat.stress(F,incomp=True))" ] }, { "cell_type": "markdown", "id": "d49731cb-b8e6-4227-9ba7-9eeeeed08a17", "metadata": {}, "source": [ "As we can see, there is a parameters `alpha_0` which is the coefficient of the spline. This can be used if the function is expressed as a linear combination of more than one splines. For example:" ] }, { "cell_type": "code", "execution_count": 17, "id": "06ae8e62-66b2-4261-bbe8-fed013b3f9a5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Material model with 2 components:\n", "Component1: splineI1\n", "Component2: splineI1\n", " ------------------------------------------------------------------\n", "Keys Value Fixed? Lower bound Upper bound \n", "------------------------------------------------------------------\n", "alpha_0 0.40 No -1.00e+01 10.00 \n", "alpha_1 0.60 No -1.00e+01 10.00 \n", "------------------------------------------------------------------\n", "\n", "[[0.6650294 0.18137165 0. ]\n", " [0.18137165 0.6801437 0. ]\n", " [0. 0. 0. ]]\n" ] } ], "source": [ "sp2 = make_interp_spline(I1,np.sin(I1-3)) #create another interpolating spline\n", "mat = pmt.MatModel('splineI1','splineI1')\n", "models = mat.models\n", "models[0].set(sp)\n", "models[1].set(sp2)\n", "params = mat.parameters\n", "params.set('alpha_0',0.4)\n", "params.set('alpha_1',0.6)\n", "mat.parameters = params\n", "print(mat, mat.parameters)\n", "print(mat.stress(F,incomp=True))" ] }, { "cell_type": "markdown", "id": "6547d340-41e1-4a47-a79f-cfa7fd16a657", "metadata": {}, "source": [ "With a spline, an interpolation is used within the range of defined values. Therefore, to track if we are extrapolating during the stress calculation, we can set the `_warn` flag to be `True`." ] }, { "cell_type": "code", "execution_count": 18, "id": "f9d65948-249f-43e5-a40e-02ed062a850e", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[0.6650294 0.18137165 0. ]\n", " [0.18137165 0.6801437 0. ]\n", " [0. 0. 0. ]]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/usr/local/lib/python3.9/site-packages/pymecht/MatModel.py:946: UserWarning: Outside the training range; be careful interpreting the results 3.8899999999999997\n", "3.0 3.1\n", " warnings.warn(w)\n" ] } ], "source": [ "models[0]._warn = True\n", "models[1]._warn = True\n", "print(mat.stress(F,incomp=True))" ] }, { "cell_type": "markdown", "id": "6102c181-5898-4f05-87b2-1b892fe2da13", "metadata": {}, "source": [ "The same procedure can be used for a bivariate spline that is a function of `I1` and `I4` defined on a rectangular grid as shown below." ] }, { "cell_type": "code", "execution_count": 19, "id": "edbfb9f5-25db-4f35-86a6-7890a4cf42fa", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Material model with 1 component:\n", "Component1: splineI1I4\n", " ------------------------------------------------------------------\n", "Keys Value Fixed? Lower bound Upper bound \n", "------------------------------------------------------------------\n", "alpha_0 1.00 No -1.00e+01 10.00 \n", "------------------------------------------------------------------\n", "\n", "[[1.27222222 0.23333333 0. ]\n", " [0.23333333 0.69166667 0. ]\n", " [0. 0. 0. ]]\n" ] } ], "source": [ "mat = pmt.MatModel('splineI1I4')\n", "I4 = np.linspace(0.9,1.1,15)\n", "I1grid, I4grid = np.meshgrid(I1,I4)\n", "Psi = (I1grid-3) + (I4grid-1)**2 + (I4grid-1)*(I1grid-3)\n", "sp = RectBivariateSpline(I1, I4, Psi.T, s=0) \n", "mat.models[0].set(sp)\n", "print(mat,mat.parameters)\n", "mat.models[0].fiber_dirs = np.array([1,0,0])\n", "print(mat.stress(F,incomp=True))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.16" } }, "nbformat": 4, "nbformat_minor": 5 }