diff --git a/examples/multi_material_tutorial/.gitignore b/examples/multi_material_tutorial/.gitignore new file mode 100644 index 00000000..87620ac7 --- /dev/null +++ b/examples/multi_material_tutorial/.gitignore @@ -0,0 +1 @@ +.ipynb_checkpoints/ diff --git a/examples/multi_material_tutorial/blocks_elements.png b/examples/multi_material_tutorial/blocks_elements.png new file mode 100644 index 00000000..f4190d0d Binary files /dev/null and b/examples/multi_material_tutorial/blocks_elements.png differ diff --git a/examples/multi_material_tutorial/boundary_conditions.jpeg b/examples/multi_material_tutorial/boundary_conditions.jpeg new file mode 100644 index 00000000..798d0190 Binary files /dev/null and b/examples/multi_material_tutorial/boundary_conditions.jpeg differ diff --git a/examples/multi_material_tutorial/optimism_tutorial.ipynb b/examples/multi_material_tutorial/optimism_tutorial.ipynb new file mode 100644 index 00000000..54af4951 --- /dev/null +++ b/examples/multi_material_tutorial/optimism_tutorial.ipynb @@ -0,0 +1,914 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b5d338a3", + "metadata": {}, + "source": [ + "# Tutorial: a multi-material simulation with optimism" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "24fd5985", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n" + ] + } + ], + "source": [ + "import numpy as onp # as in \"old\" numpy, to distinguish it from Jax's numpy.\n", + "\n", + "from optimism.JaxConfig import *\n", + "from optimism import EquationSolver\n", + "from optimism import FunctionSpace\n", + "from optimism import Interpolants\n", + "from optimism.material import J2Plastic\n", + "from optimism.material import Neohookean\n", + "from optimism import Mechanics\n", + "from optimism import Mesh\n", + "from optimism import Objective\n", + "from optimism import QuadratureRule\n", + "from optimism import SparseMatrixAssembler\n", + "from optimism import VTKWriter" + ] + }, + { + "cell_type": "markdown", + "id": "9c7de108", + "metadata": {}, + "source": [ + "## Overview of the problem\n", + "\n", + "In this problem, we'll subject a rectangular slab to plane strain uniaxial tension. As shown in the figure, the slab is made of two equal-size layers made of different materials. One material will use a hyperelastic material model, while the other will use a $J_2$ plasticity model." + ] + }, + { + "cell_type": "markdown", + "id": "37c03013", + "metadata": {}, + "source": [ + "*Figure here*" + ] + }, + { + "cell_type": "markdown", + "id": "b1c62590", + "metadata": {}, + "source": [ + "## Set up the mesh" + ] + }, + { + "cell_type": "markdown", + "id": "a8b6f7ba", + "metadata": {}, + "source": [ + "As in any finite element code, first we need to set up the geometry of the problem with a mesh. Optimism can read in meshes from the Exodus II format, but it also provides utilities for generating structured meshes for simple problems like this one. Let's generate a rectangular mesh with width `w` and length `L`. We'll also specify that it has 5 nodes in the $x$-direction (hence 4 elements) and 15 nodes in the $y$-direction." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "b7750f13", + "metadata": {}, + "outputs": [], + "source": [ + "# set up the mesh\n", + "w = 0.1\n", + "L = 1.0\n", + "nodesInX = 5 # must be odd\n", + "nodesInY = 15\n", + "xRange = [0.0, w]\n", + "yRange = [0.0, L]\n", + "mesh = Mesh.construct_structured_mesh(nodesInX, nodesInY, xRange, yRange)" + ] + }, + { + "cell_type": "markdown", + "id": "f18bd7ea", + "metadata": {}, + "source": [ + "The mesh object has an attribute called `blocks` which lets us group elements together to specify different materials on them. Let's see what our mesh has:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "bb1bd192", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'block_0': DeviceArray([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,\n", + " 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,\n", + " 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,\n", + " 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47,\n", + " 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59,\n", + " 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71,\n", + " 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83,\n", + " 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95,\n", + " 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107,\n", + " 108, 109, 110, 111], dtype=int64)}\n" + ] + } + ], + "source": [ + "print(mesh.blocks)" + ] + }, + { + "cell_type": "markdown", + "id": "2bea405c", + "metadata": {}, + "source": [ + "`blocks` is a standard Python dictionary object (i.e. a `dict`), where the keys are strings that lets us give meaningful names to the element blocks. The values are Jax-numpy arrays that contain the indices of the elements in the block. The `construct_structured_mesh` function returns a mesh with just one block. We want two equal-sized blocks for our problem like in the figure, so let's modify the mesh.\n", + "\n", + "First, we'll define a simple helper function that takes in the vertices of a triangular element and returns the coordinates of the centroid." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "ec126372", + "metadata": {}, + "outputs": [], + "source": [ + "def triangle_centroid(vertices):\n", + " return np.average(vertices, axis=0)" + ] + }, + { + "cell_type": "markdown", + "id": "d88821e7", + "metadata": {}, + "source": [ + "We'll loop over the element node connectivity table, which is the `mesh.conns` object, extract the coordinates of its vertices, and use the `triangle_centroid` object on them to determine if the element is in the left or right layer. We will store the results in a list called `blockIds`:" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "a35d7cdd", + "metadata": {}, + "outputs": [], + "source": [ + "# initialize the block IDs of all elements to a dummy value (-1)\n", + "blockIds = -1*onp.ones(Mesh.num_elements(mesh), dtype=onp.int64)\n", + "for e, t in enumerate(mesh.conns):\n", + " vertices = mesh.coords[t,:]\n", + " if triangle_centroid(vertices)[0] < w/2:\n", + " blockIds[e] = 0\n", + " else:\n", + " blockIds[e] = 1\n", + "# check that every element has gotten an ID\n", + "assert(onp.all(blockIds != -1))" + ] + }, + { + "cell_type": "markdown", + "id": "b87eed4b", + "metadata": {}, + "source": [ + "This will mark an element as block 0 if it's centroid is on the left hand side of the slab, and as block 1 if it's on the other side. Now, let's make the `dict` that we want to attach to the mesh object." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "e245933d", + "metadata": {}, + "outputs": [], + "source": [ + "blocks = {'left layer': np.flatnonzero(np.array(blockIds) == 0),\n", + " 'right layer': np.flatnonzero(np.array(blockIds) == 1)}" + ] + }, + { + "cell_type": "markdown", + "id": "61954ea8", + "metadata": {}, + "source": [ + "Now we can make use of a function that takes in the original mesh (with one block) and the block IDS we just created, and returns a new mesh that is the same as the old one except that the blocks have been updated." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "e2a1c548", + "metadata": {}, + "outputs": [], + "source": [ + "mesh = Mesh.mesh_with_blocks(mesh, blocks)" + ] + }, + { + "cell_type": "markdown", + "id": "bcb47738", + "metadata": {}, + "source": [ + "Let's check to make sure this process worked. To do this, we'll make use of optimism's output facilities. Optimism provides a class called `VTKWriter`, that, as the name suggests, writes data to VTK files which can be visualized in ParaView (and several other visualization tools). First, we instantiate a VTKWriter object, giving it the base of the filename (the name that will be suffixed with the \".vtk\" extension)." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "ea13862b", + "metadata": {}, + "outputs": [], + "source": [ + "writer = VTKWriter.VTKWriter(mesh, baseFileName='check_problem_setup')" + ] + }, + { + "cell_type": "markdown", + "id": "619fe562", + "metadata": {}, + "source": [ + "Next, we can start adding fields. In this case, we only have one field - the block ID numbers - which is a scalar field of integer type. Finally, once our data is loaded into the writer, we call the `write()` method to tell it to write the VTK file to disk." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "42104d3f", + "metadata": {}, + "outputs": [], + "source": [ + "writer.add_cell_field(name='block_id', cellData=blockIds,\n", + " fieldType=VTKWriter.VTKFieldType.SCALARS,\n", + " dataType=VTKWriter.VTKDataType.INT)\n", + "writer.write()" + ] + }, + { + "cell_type": "markdown", + "id": "607566b1", + "metadata": {}, + "source": [ + "This is what we get when we open Paraview and visualize the `block_id` field. If you try it, you should see something similar." + ] + }, + { + "cell_type": "markdown", + "id": "50aa988b", + "metadata": {}, + "source": [ + "![paraview shwoing blocks in mesh](blocks_elements.png)" + ] + }, + { + "cell_type": "markdown", + "id": "f0637278", + "metadata": {}, + "source": [ + "Success! This is what were shooting for." + ] + }, + { + "cell_type": "markdown", + "id": "bddeee2a", + "metadata": {}, + "source": [ + "## Essential boundary conditions\n", + "We're going to make one more modification to the mesh. Looking again at the problem setup figure, we can see that we need to apply boundary conditions on the bottom, left, and top boundary edges of the slab. Similar to the `blocks` attribute, `mesh` has a `nodeSets` dictionary that maps names to sets of nodes. We'll make the nodesets we need by performing range queries over the nodal coordinates:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "4d67689c", + "metadata": {}, + "outputs": [], + "source": [ + "nodeTol = 1e-8\n", + "nodeSets = {'left': np.flatnonzero(mesh.coords[:,0] < xRange[0] + nodeTol),\n", + " 'right': np.flatnonzero(mesh.coords[:,0] > xRange[1] - nodeTol),\n", + " 'bottom': np.flatnonzero(mesh.coords[:,1] < yRange[0] + nodeTol),\n", + " 'top': np.flatnonzero(mesh.coords[:,1] > yRange[1] - nodeTol)}\n", + "\n", + "# create a copy of the mesh that has the nodeSets in it\n", + "mesh = Mesh.mesh_with_nodesets(mesh, nodeSets)" + ] + }, + { + "cell_type": "markdown", + "id": "413f82f9", + "metadata": {}, + "source": [ + "Now we're going to register the essential boundary conditions so that the optimism solvers will know how to deal with them. This is done with an `EssentialBC` object. Each `EssentialBC` represents a boundary condition on one field component of a node set. As en example, let's create one to represent the $x$-component displacement boundary condition on the nodes of the left edge:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "ebedd68d", + "metadata": {}, + "outputs": [], + "source": [ + "ebcLeft = Mesh.EssentialBC(nodeSet='left', field=0)" + ] + }, + { + "cell_type": "markdown", + "id": "3ffe98a4", + "metadata": {}, + "source": [ + "This is one boundary condition; we have three essential boundary conditions in total to apply. The thing to do is to collect them all in a list. So the code looks like it would in a real application, we'll ignore the `ebcLeft` we just created and create all of the essential boundary conditions in one statement:" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "70344154", + "metadata": {}, + "outputs": [], + "source": [ + "EBCs = [Mesh.EssentialBC(nodeSet='left', field=0),\n", + " Mesh.EssentialBC(nodeSet='bottom', field=1),\n", + " Mesh.EssentialBC(nodeSet='top', field = 1)]" + ] + }, + { + "cell_type": "markdown", + "id": "9f671f9e", + "metadata": {}, + "source": [ + "Next, we create a `DofManager` object. What is this for? It's a class to help us index into our nodal arrays, keeping track of which degrees of freedom have essential boundary conditions. It's purpose will become clearer when we use it later." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "a34991d4", + "metadata": {}, + "outputs": [], + "source": [ + "fieldShape = mesh.coords.shape\n", + "dofManager = Mesh.DofManager(mesh, fieldShape, EBCs)" + ] + }, + { + "cell_type": "markdown", + "id": "0c4b89c1", + "metadata": {}, + "source": [ + "The variable `fieldShape` tells `dofManager` what the array of the unknowns will look like. In solid mechanics, the unknown field is the displacement, which happen to have the same shape as the nodal coordinates, so the `mesh.coords` array is a convenient place to grab this information. You could also manually enter `(nNodes, 2)`, where `nNodes` is the total number of nodes in the mesh, and 2 is the number of components of displacement in plane strain.\n", + "\n", + "We use the vis tools again to check that we've done the essential boundary condition specification correctly. First, we'll use a little trick to turn the boundary conditions into a viewable nodal field. We'll use an attribute of the `DofManager` class called `isBc`. This is just an array of size `fieldShape` that indicates whether each dof has an essential boundary condition. Numpy can cast this to an array of integers (more precisely, the `int64` type in numpy) with values 0 or 1 which can be plotted in Paraview. The `dataType` argument is different now; for block ID, it was a scalar field, but for boundary conditions, we want it to be a vector field (one component for each displacement component)." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "3d92324c", + "metadata": {}, + "outputs": [], + "source": [ + "bcs = np.array(dofManager.isBc, dtype=np.int64)\n", + "writer.add_nodal_field(name='bcs', nodalData=bcs, fieldType=VTKWriter.VTKFieldType.VECTORS,\n", + " dataType=VTKWriter.VTKDataType.INT)\n", + "writer.write()" + ] + }, + { + "cell_type": "markdown", + "id": "1103ae8f", + "metadata": {}, + "source": [ + "Note that `writer` still refers to the same `VTKWriter` object as before. A VTKWriter object is always associated with the same filename, so when we add a new field and then call the `write()` method, it will overwrite the previous VTK file. Indeed, if you open `check_problem_steup.vtk`, you'll see that it now contains two output fields, \"block_id\" and \"bcs\".\n", + "\n", + "Contour plots of the $x$- and $y$-components of the \"bcs\" field are shown below:" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "ceeb6e96", + "metadata": {}, + "source": [ + "![contour plots of boundary condition fields](boundary_conditions.jpeg)" + ] + }, + { + "cell_type": "markdown", + "id": "f4ee2543", + "metadata": {}, + "source": [ + "The first plot shows all nodes with boundary conditions on the $x$ component of the displacement. We see that the left edge has a value of 1 (which means that the boundary condition flag is \"True\" there), and every other node has a value of 0, which means they are unrestrained. This is exactly what we want. The $y$ component plot also confirms that the top and bottom edges correctly have their boundary conditions. Of course, the top edge has an *inhomogeneous* boundary condition. We'll enforce the actual value of this boundary condition later." + ] + }, + { + "cell_type": "markdown", + "id": "f8cef770", + "metadata": {}, + "source": [ + "## Build the function space\n", + "The next step is to create a representation of the discrete function space to help us do things like interpolate our fields and do calculus on them. The first ingredient we need is a quadrature rule. In optimism, quadrature rules are specified by the highest degree polynomial that they can exactly integrate. A smart way to do this is to choose it in relation to $p$, the polynomial degree of our interpolation." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "68a56fe0", + "metadata": {}, + "outputs": [], + "source": [ + "p = mesh.masterElement.degree" + ] + }, + { + "cell_type": "markdown", + "id": "d0b6b351", + "metadata": {}, + "source": [ + "In the linearized theory of solid mechanics, the natural trial/test function space is $H^1$, because the operator contains products of gradients of the displacement. Since our interpolation is of degree $p$, the displacement gradient is of degree $p-1$, and the inner product of gradients is of degree $2(p-1)$. Thus, we choose this as our quadrature rule precision to avoid under-integrating our operator:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "cf915ada", + "metadata": {}, + "outputs": [], + "source": [ + "# create the function space\n", + "quadRule = QuadratureRule.create_quadrature_rule_on_triangle(degree=2*(p - 1))" + ] + }, + { + "cell_type": "markdown", + "id": "fbb9d8e4", + "metadata": {}, + "source": [ + "The benefit of specifying our quadrature rule in this way is that if we later decide to modify the mesh to use higher-order elements, the quadrature rule will be updated automatically. This helps keep us out of trouble with hourglassing problems. Note that our operator is *nonlinear*, so the quadrature rule won't integrate it exactly, but the accuracy requirement of the linearized theory is a good rule of thumb.\n", + "\n", + "With the quadrature rule (and the mesh), we can now construct the function space:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "6c36cd8b", + "metadata": {}, + "outputs": [], + "source": [ + "fs = FunctionSpace.construct_function_space(mesh, quadRule)" + ] + }, + { + "cell_type": "markdown", + "id": "f79a4bc8", + "metadata": {}, + "source": [ + "We'll see this object in operation later when we set up our energy function." + ] + }, + { + "cell_type": "markdown", + "id": "d44fd1a0", + "metadata": {}, + "source": [ + "## Material models\n", + "Next, we instantiate the material models for the slab. For the left side, we'll choose a Neo-Hookean material. Material models in optimism take their material parameters in the form of a dictionary. For Neo-Hookean, the required parameters are the elastic modulus and the Poisson's ratio (both taken about the undeformed state)." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "2c1f89f0", + "metadata": {}, + "outputs": [], + "source": [ + "# create the material model for the left side\n", + "props = {'elastic modulus': 5.0,\n", + " 'poisson ratio': 0.25}\n", + "leftMaterialModel = Neohookean.create_material_model_functions(props)" + ] + }, + { + "cell_type": "markdown", + "id": "66eab9e6", + "metadata": {}, + "source": [ + "TODO: The property names are not documented yet. For now, you can find them by inspecting the code in optimism/materials." + ] + }, + { + "cell_type": "markdown", + "id": "176065d4", + "metadata": {}, + "source": [ + "Now we'll instantiate the other material model for the right-hand side. We'll pick a $J_2$ plasticity model, which is a bit more interesting (and thus has more parameters)." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "24778eaf", + "metadata": {}, + "outputs": [], + "source": [ + "E = 10.0\n", + "nu = 0.25\n", + "Y0 = 0.01*E\n", + "H = E/100\n", + "props = {'elastic modulus': E,\n", + " 'poisson ratio': nu,\n", + " 'yield strength': Y0,\n", + " 'hardening model': 'linear',\n", + " 'hardening modulus': H}\n", + "rightMaterialModel = J2Plastic.create_material_model_functions(props)\n" + ] + }, + { + "cell_type": "markdown", + "id": "19011b92", + "metadata": {}, + "source": [ + "The meaning of the parameters is clear from the keys. There are several hardening models currently available, such as linear hardening, a version of power law hardening, and the Voce exponential saturation law. We'll keep it simple and just use linear hardening.\n", + "\n", + "For multi-material simulations, we must create a dictionary of each material that maps the name of each element block to the material model object for that block, like this:" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "aa2b0a4e", + "metadata": {}, + "outputs": [], + "source": [ + "materialModels = {'left layer': leftMaterialModel, 'right layer': rightMaterialModel}\n" + ] + }, + { + "cell_type": "markdown", + "id": "3f1723e8", + "metadata": {}, + "source": [ + "## Write the energy function to minimize" + ] + }, + { + "cell_type": "markdown", + "id": "94039b0f", + "metadata": {}, + "source": [ + "Numerical solutions to PDEs are obtained in optimism by minimizing an objective function, which may be thought of intuitively as en energy. In fact, for hyperelastic materials, the objective function *is* the total energy. For history-dependent materials, one can often write an incremental functional that is minimized by the displacement field that carries the body from one discrete time step to the next. We're going to write the function for our problem now.\n", + "\n", + "There is a tool in the `Mechanics` module that will do most of the work for us. Let's call it first and explain its output afterwards." + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "e20f5769", + "metadata": {}, + "outputs": [], + "source": [ + " \n", + "# mechanics functions\n", + "mechanicsFunctions = Mechanics.create_multi_block_mechanics_functions(\n", + " fs, mode2D=\"plane strain\", materialModels=materialModels)\n" + ] + }, + { + "cell_type": "markdown", + "id": "35578c4a", + "metadata": {}, + "source": [ + "The `create_multi_block_mechanics_functions` writes some functions for us to help write the problem in the right form. The most important part of the output is `mechanicsFunctions.compute_strain_energy(...)`, which writes all the loops over the blocks, elements, and quadrature points that you would need to compute the energy from the nodal displacements. (Now we finally see the `FunctionSpace` object `fs` in action).\n", + "To use this new function, we'll invoke it like this:\n", + "```python\n", + "mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", + "```\n", + "where `U` is the array of the degrees of freedom (the nodal displacements), and `internalVariables` is an array of the internal variables at every quadrature point. It has shape `(ne, neqp, nint)`, with `ne` being the number of elements, `neqp` the number of quadrature points per element, and `nint` the number of internal variables for the material model. (NOTE: For multi-material simulations, this is currently sized such that `nint` is the *maximum* number of internal variables over all materials, so that the array is padded for materials using fewer internal variables. We may remove this restriction in the future to improve memory efficiency).\n", + "\n", + "All of the minimization solvers in optimism require the objective function to have a signature like this:\n", + "```python\n", + "energy_function(Uu, p)\n", + "```\n", + "where `Uu` is the array of unknowns, and `p` is called the parameter set. The parameter set essentially holds all of the information that is needed to specify the problem but *isn't* the set of unknowns. These are things like values of the boundary conditions and the internal variables, as well as some other categories. Parameter sets are constructed by calling the `Params` function in the `Objective` module. This is to help organize them in certain categories that the solver needs to be aware of, such as which ones have derivatives taken with respect to inside the solver. We need to cast our energy function in this form. Let's write it like this and work out what the intervening steps must be:\n", + "\n", + "```python\n", + "def energy_function(Uu, p):\n", + " U = create_field(Uu, p)\n", + " internalVariables = p.state_data\n", + " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", + "```\n", + "On the first line, we use the parameter set to pass from the array of unknowns to the full array of nodal displacements. This means we need to fill in the values of the inhomogeneous boundary conditions. Next, we pull out the internal variables from the parameter set. Finally, we use the canned `compute_strain_energy(...)` function with these variables in order to compute the total energy.\n", + "\n", + "The inhomogeneous boundary condition part is handled like so:" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "217248f1", + "metadata": {}, + "outputs": [], + "source": [ + "# helper function to fill in nodal values of essential boundary conditions\n", + "def get_ubcs(p):\n", + " appliedDisp = p[0]\n", + " EbcIndex = (mesh.nodeSets['top'], 1)\n", + " V = np.zeros(fieldShape).at[EbcIndex].set(appliedDisp)\n", + " return dofManager.get_bc_values(V)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "c578d78c", + "metadata": {}, + "source": [ + "We will store the applied displacement in the first slot of the parameter set. In line 4 above we extract it. Then we make an array of the same size as the nodal displacements, set it to zero, and replace the values in the DOFs on the top edge with the value of the applied displacement.\n", + "\n", + "Now we can write the `create_field` function shown above in the proposed `energy_function`:" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "id": "2755c6db", + "metadata": {}, + "outputs": [], + "source": [ + "# helper function to go from unknowns to full DoF array\n", + "def create_field(Uu, p):\n", + " return dofManager.create_field(Uu, get_ubcs(p))\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "id": "34511c5d", + "metadata": {}, + "source": [ + "The pieces are finally in place to define the energy function for our problem:" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "be482478", + "metadata": {}, + "outputs": [], + "source": [ + "# write the energy to minimize\n", + "def energy_function(Uu, p):\n", + " U = create_field(Uu, p)\n", + " internalVariables = p.state_data\n", + " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "5fab41f4", + "metadata": {}, + "source": [ + "## Set up the optimization solver\n", + "We have an objective function - `energy_function`, which we will hand to a routine that will find the unknowns displacements that minimize it. In this section, we specficy which optimization solver we want to use, and tell it how to work. \n", + "\n", + "We will use the Steighaug trust region method. This method uses linear conjugate gradient iterations as part of its algorithm, which in turn need to be preconditioned in order to be effective. Currently, the only available preconditioner in optimism is a Cholesky factorization on the stiffness matrix. We need to intruct the solver how to assemble the stiffness matrix like this:" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "id": "cd2bb5dd", + "metadata": {}, + "outputs": [], + "source": [ + "# Tell objective how to assemble preconditioner matrix\n", + "def assemble_sparse_preconditioner_matrix(Uu, p):\n", + " U = create_field(Uu, p)\n", + " internalVariables = p.state_data\n", + " elementStiffnesses = mechanicsFunctions.compute_element_stiffnesses(U, internalVariables)\n", + " return SparseMatrixAssembler.assemble_sparse_stiffness_matrix(\n", + " elementStiffnesses, mesh.conns, dofManager)\n" + ] + }, + { + "cell_type": "markdown", + "id": "2b5da4b6", + "metadata": {}, + "source": [ + "We see once again that `mechanicsFunctions` provides a helper function. In this case, the function `compute_element_stiffnesses` takes the same inputs as the energy function, but instead of returning the total energy, it returns an array containing the stiffness matrix for each element. The elemental stiffness matrices are the Hessians of the total energy in each element, and automatic differentiation is again used to perform this calculation. The `assemble_sparse_stiffness_matrix(...)` function takes these elemental stiffness matrices and contructs the system's stiffness matrix using a sparse matrix data structure (from SciPy). We tell the solver how to use this capability by building something called a `PrecondStrategy`: " + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "id": "b1bc05ed", + "metadata": {}, + "outputs": [], + "source": [ + "precondStrategy = Objective.PrecondStrategy(assemble_sparse_preconditioner_matrix)\n", + "\n", + "# solver settings\n", + "solverSettings = EquationSolver.get_settings(use_preconditioned_inner_product_for_cg=True)\n" + ] + }, + { + "cell_type": "markdown", + "id": "d4b55505", + "metadata": {}, + "source": [ + "Finally, we can specify custom settings for the solver is we wish." + ] + }, + { + "cell_type": "markdown", + "id": "3e10766a", + "metadata": {}, + "source": [ + "## Solve!" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "28b893dd", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Assembling preconditioner 0\n", + "min, max diagonal stiffness = 0.9999999999999998 1.0000000000000002\n", + "Factorizing preconditioner\n", + "num warm start cg iters = 1\n", + "Assembling preconditioner 0\n", + "min, max diagonal stiffness = 0.9875365254464693 1.0110599625088337\n", + "Factorizing preconditioner\n", + "\n", + "Initial objective, residual = 0.04170630570636782 1.5580526760827706e-05\n", + "obj= -1.47133479e-10 , model obj= -1.47130159e-10 , res= 5.10352716e-10 , model res= 7.43900588e-21 , cg iters= 1 , tr size= 2.0 , interior , accepted= True\n", + "\n" + ] + } + ], + "source": [ + "# initialize unknown displacements to zero\n", + "Uu = dofManager.get_unknown_values(np.zeros(fieldShape))\n", + "\n", + "# set initial values of parameters\n", + "appliedDisp = 0.0\n", + "state = mechanicsFunctions.compute_initial_state()\n", + "p = Objective.Params(appliedDisp, state)\n", + " \n", + "# Construct an objective object for the equation solver to work on\n", + "objective = Objective.ScaledObjective(energy_function, Uu, p, precondStrategy=precondStrategy)\n", + " \n", + "# increment the applied displacement\n", + "appliedDisp = L*0.01\n", + "p = Objective.param_index_update(p, 0, appliedDisp)\n", + "\n", + "# Find unknown displacements by minimizing the objective\n", + "Uu = EquationSolver.nonlinear_equation_solve(objective, Uu, p, solverSettings)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "2b1968f3", + "metadata": {}, + "source": [ + "## Post-process\n", + "The last step is to write out the simulation results so that we can use them. We've already seen a few examples of generating VTK output. Let's create another `VTKWriter` for the simulation results. Adding the displacement field is straightforward:" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "e14a722a", + "metadata": {}, + "outputs": [], + "source": [ + "# write solution data to VTK file for post-processing\n", + "writer = VTKWriter.VTKWriter(mesh, baseFileName='uniaxial_two_material')\n", + "\n", + "U = create_field(Uu, p)\n", + "writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "da5415dd", + "metadata": {}, + "source": [ + "Let's also show an example of plotting quadrature field data. A commonly needed output is the stress field. We make use of another function in `mechanicsFunctions` to help us. However, before we write out any quadrature field data, we should update the internal variables. Currently, `state` still holds the initial conditions of the internal variables. That is, the solver finds the equilibrium displacement field, but doesn't change `state` in place. This is again due to Jax's functional programming paradigm. The following function call returns the updated internal variables using the new displacement field and the old internal variables as inputs:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e46298c0", + "metadata": {}, + "outputs": [], + "source": [ + "# update the state variables in the new equilibrium configuration\n", + "state = mechanicsFunctions.compute_updated_internal_variables(U, p.state_data)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "d23fc4a2", + "metadata": {}, + "source": [ + "We make use of another of the `mechanicsFunctions` member functions to get the stress field, using the updated internal variables. Note that we don't pay the cost of iteratively solving the history-dependent material model response again. This function expects that the internal variables are already updated when it is called, and simply uses theses values in the evaluation of the material model energy density functions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3b3e743", + "metadata": {}, + "outputs": [], + "source": [ + "energyDensities, stresses = mechanicsFunctions.compute_output_energy_densities_and_stresses(\n", + " U, state)\n" + ] + }, + { + "cell_type": "markdown", + "id": "b4221fcc", + "metadata": {}, + "source": [ + "As the left-hand side of this assignment implies, the scalar energy density field (at every quadrature point) is returned in addition to the stresses. In fact, only the scalar-valued energy density function for each material is implemented in each material model. The stress field is the derivative of the energy density function with respect to the displacement gradient, and optimism uses Jax's automatic differentiation to compute this derivative. When computing this derivative, evaluating the energy density at the same time is essentially free, so the above function coputes them both.\n", + "\n", + "The `stresses` array contains the stress tensor for every quadrature point of every element. There is no interpolation associated with a quadrature field, so in order to visualize the stresses, this must be remedied, either on the physics solver side or in the visualization software. One simple way to accomplish this is to compute element averages of the field and plot it as a cell-based field (instead of a nodal field). Optimism provides a method to do this:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0d3a7bb", + "metadata": {}, + "outputs": [], + "source": [ + "cellStresses = FunctionSpace.project_quadrature_field_to_element_field(fs, stresses)\n", + "writer.add_cell_field(name='stress', cellData=cellStresses,\n", + " fieldType=VTKWriter.VTKFieldType.TENSORS)\n", + "\n", + "writer.write()\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "2558d9ca", + "metadata": {}, + "source": [ + "The above code can be wrapped in a loop that will compute the response at multiple time steps. Note that before moving to the next time step, you should update the internal variables in the parameter set:" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "7e0acf41", + "metadata": {}, + "outputs": [], + "source": [ + "# update the state variables in the new equilibrium configuration\n", + "p = Objective.param_index_update(p, 1, state)" + ] + } + ], + "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.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}