Faim Finite Element Analysis Manual

The solvers and n88modelgenerator are the non-open sourced components of Faim and require a license file to run. They are bundled together into installers that can be downloaded from https://numerics88.com/downloads/ .

The Linux installer will be named something like faim-9.0-linux.sh. It can be installed with the following command for a single user:

For a system-wide installation, add sudo:

If you want to specify a particular install location, use the --prefix flag, for example

The Windows installer will be named something like faim-9.0-windows.exe. Simply run it to install.

The Mac installer is a disk image with a name like faim-9.0-windows.dmg. Double-click to mount it, and then drag the Faim icon to the Applications folder, or to any other convenient location.

For all operating systems, it is possible to install Faim system-wide using an administrator account or for a single user, for which administrator rights are not needed. It is also possible to install it multiple times on a single machine, so that each user has their own installation. The installation directory can also be freely renamed and moved around.

To run the solvers and n88modelgenerator, you may use any of the following methods:

To receive a license file, you must run one of the Faim programs with the --license_check option, and send the resulting UUID to skboyd@ucalgary.ca . For example,

When you receive a license file, which is a short text file, you can copy it to one of several possible locations where the software can find it. These locations are:

You may have multiple license files installed at once.

The solvers and n88modelgenerator are most conveniently used together with a collection a command-line utilities called n88tools. n88tools is implemented in Python, and depends on the library vtkbone. Additionally, vtkbone can be used to create custom model generation and processing scripts in Python, as will be discussed in a subsequent chapter.

n88tools and vtkbone are open-source. If you wish, you can download the source code from https://github.com/Numerics88/ and compile them yourself, using whatever installation of Python is most convenient for you. Most users however, will not want to go through the hassle of compiling them. The quickest way to install them is via Anaconda Python, which is supported on Linux, Windows, and macOS. An advantage of using Anaconda Python is that the very large repository of mathematical and scientific Python packages available in Anaconda Python can be used together with Faim. To install n88tools and vtkbone in Anaconda Python, follow these steps:

To use n88tools, you simply activate the Anaconda environment. On Linux and macOS, this is done with

On Windows, the command is

You will have to activate the environment in each Terminal (or Command Prompt) in which you want to use Faim.

ParaView can be used for interactive rendering of Faim finite elements models. For this purpose, Numerics88 provides plugins for ParaView. These allow ParaView to open n88model files (as well as Scanco AIM files). The plugins can be downloaded from http://numerics88.com/downloads/ . You can unzip the plugins into any convenient directory, and you can move and rename the plugin directory as you like.

To load the plugins, open ParaView, and select the menu item Tools › Manage Plugins. The Plugin Manager will appear, as shown in Figure 1.

Click Load New…​ and navigate to the location of the plugins. Select one of the plugins, and then click OK. You will need to repeat this for each of the plugins, in particular for

They are now loaded, and you should be able to read to corresponding file formats. You may want to select Auto Load for each plugin so that you don’t need to load them every time you open Paraview. To do this, click the drop-down arrow beside the plugin in the Plugin Manager, and make sure that Auto Load is selected. This is shown in Figure 2.

A PDF version of this manual can be downloaded from http://numerics88.com/documentation/ .

Example data files, together with scripts for the tutorials in this manual, can be downloaded from http://numerics88.com/downloads/ .

The first step is to segment your 3D medical image data. Although in principle any type of 3D image data can be analyzed, we focus on the use of micro-CT images of bone here. During segmentation, bone and other relevant structures in the CT image are labelled and identified. Segmentation is a very broad and complex field. Because of this, Numerics88 software does not attempt to provide all solutions. Instead, we recommend that you use the segmentation tool provided by the manufacturer of your medical image equipment. This has the advantage that it will be well-integrated into the scanning workflow. In the case of Scanco systems, the tool is Image Processing Language (IPL); for Skyscan systems, CT-analyzer provides segmentation functionality.

In cases where your segmentation problem is unusually tricky and beyond the abilities of the manufacturer’s tools, there are many third party solutions available.

Whichever tool you use to perform the segmentation, the result should be a 3D image file with integer values. These integer values are referred to as "Material IDs", because they will subsequently be associated with abstract mathematical material models. Material ID zero indicates background material (often air, but sometimes also other material, such as marrow in bone) that can be assumed to have negligible stiffness. Voxels with material ID zero will not be converted to elements in the FE model. Other values can be associated freely with particular materials or material properties as required.

Model generation is the process of converting a segmented image to a geometrical representation suitable for finite element analysis, complete with suitable constraints such as displacement boundary conditions and applied loads.

A number of steps are required for complete model generation.

Faim provides two tools for model generation:

The finite element solver takes the input model and calculates a solution of the displacements and forces on all the nodes, subject to the specified constraints.

As finite element models derived from micro-CT can be very large, the memory efficiency and speed of the solver are important. Faim uses a mesh-free preconditioned conjugate gradient iterative solver. This type of solver is highly memory efficient.

Post-processing and visualization can take many different forms, depending on the goals of your analysis. For example, you may want to identify locations of stress variations, in which case, a visualization tool is likely the best approach. Or, you may want to obtain a basic output such as an overall stiffness. A number of common post-processing operations are performed by Faim.

Visualization is important even in cases where ultimately you want one or more quantitative values. Visualization provides a conceptual overview of the solution, and may allow one to quickly identify incorrect or invalid results due to some deficiency in the input model. Many visualization software packages exist. If you don’t have an existing favourite package, we recommend the use of ParaView, which is an open-source solution from Kitware. The Faim distribution includes plug-ins for ParaView which allow ParaView to directly open Faim file types.

A uniaxial test is a compression (or tension) test with a fixed displacement or force applied along the z axis. Two boundary conditions are applied:

A uniaxial test is shown in Figure 5, along with other types of compression tests. The coordinate system shown in the figure is the Test Frame. A uniaxial test is distinguished from an axial test (described below) by the fact that nodes on the top and bottom surfaces are unconstrained laterally. This corresponds to zero contact friction.

The amount of displacement applied to the top surface can be specified by either the normal_strain or displacement parameters.

An axial test is similar to a uniaxial test, with the additional constraint that nodes on the top and bottom surfaces are also laterally fixed (i.e. no movement of these nodes is permitted in the x-y directions). This corresponds to 100% contact friction.

An axial test is shown in Figure 5.

A confined test is similar to an axial test, with the addition that nodes on the side surfaces of the image volume are constrained laterally (i.e. no movement of these nodes is permitted in the x-y directions).

A confined test is shown in Figure 5.

In a symmetric shear test, the side surfaces are angularly displaced to correspond to a given shear strain. This is shown in Figure 6. The only parameter of relevance to a symmetric shear test is shear_strain.

In a directional shear test, the top surface is displaced laterally, as shown in Figure 7. The direction and degree of displacement are set by the parameter shear_vector, which gives the x,y shear displacement in the Test Frame. Typically shear_vector is unitless, and the actual displacement is the shear vector scaled by the vertical height of the model (i.e. the image extent in the z direction in the Test Frame). If however scale_shear_to_height is set to off, then shear_vector is taken to have absolute units of length.

In a bending test, the top and bottom surfaces are rotated in opposite directions, as shown in Figure 8. The rotation axes are defined by a neutral axis, given by a point in the x,y plane (shown as C) and an angle in the x,y plane. (The coordinate system shown in the figure is the Test Frame.) The neutral axis, projected onto the top and bottom surfaces, is indicated in blue in Figure 8. The point and angle defining the neutral axis can be set in n88modelgenerator with the parameters central_axis and neutral_axis_angle respectively. The central axis (point C) can be specified as a numerical x,y pair, or as the center of mass or the center of the data bounds; center of mass is the default setting. The default angle of the neutral axis is 90º, parallel to the y-axis. The amount of rotation can be specified with the parameter bending_angle. The top and the bottom surfaces are each rotated by half of this value. Positive rotation is defined as in Figure 8.

In a bending test, the top and bottom surfaces are laterally constrained; no x,y motion of the nodes on these surfaces is permitted.

In a torsion test, the top surface is rotated about the center axis, as shown in Figure 9. The angle of rotation is given by twist_angle. In the figure positive rotation is shown. The central axis lies parallel to the z axis in the test frame; its position can be specified with the parameter central_axis. The default position of the central axis is passing through the center of mass.

Linear elastic materials obey Hooke’s Law,

where ε_xx,ε_yy,ε_zz are the engineering normal strains along the axis directions, γ_yz,γ_zx,γ_xy are the engineering shear strains, and σ_ij are the corresponding stresses. S is the compliance matrix. The inverse of S gives σ in terms of ε, and is the stress-strain matrix C, also sometimes called the stiffness matrix.

An isotropic material has the same mechanical properties in every direction. Elastic isotropic materials are characterized by only two independent parameters: Young’s modulus (E) and Poisson’s ratio (ν). The compliance matrix is

For n88modelgenerator the isotropic elastic parameters are set with youngs_modulus and poissons_ratio. The default values are 6829 MPa and 0.3, as reported by MacNeil and Boyd (2008).

An isotropic material definition is the default, and will be assumed if no other material type is specified.

Orthotropic elastic materials have two or three mutually orthogonal twofold axes of symmetry. The compliance matrix is

where

Note that the ν_ij are not all independent. The compliance matrix must be symmetric, hence

Therefore only 9 independent parameters are required to fully specify an orthotropic material.

For n88modelgenerator the orthotropic elastic parameters are set with orthotropic_parameters. Whenever this parameter is specified, an orthotropic material definition will be used.

The most general form of an elastic material is anisotropic, for which we allow the stress-strain matrix, or the compliance matrix, to have arbitrary values. The stress-strain matrix is symmetric, therefore an anisotropic material requires 21 parameters to define. To define an anisotropic material with n88modelgenerator, it is necessary to use a material definitions file. See below.

Elastoplastic material behaviour is an idealisation of a particular kind of nonlinear stress-strain relationship in which we divide the stress-strain relationship into two clearly distinguished regions: (1) an elastic region, which is identical to the linear elastic material behaviour described above, and, (2) a purely plastic region. For a material subject to a uniaxial deformation, this is sketched in Figure 10. In Figure 10A, we see the initial elastic region, which in this sketch is tensile, with slope equal to the Young’s modulus. This region continues until the stress reaches the yield strength, Y. For further increases in strain beyond this point, the stress ceases to increase, and remains constant at the yield strength. The incremental strain in this region is referred to as the plastic strain. Although we must perform work to produce a plastic strain, the plastic strain is irreversible, as can be seen in Figure 10B. In B, after applying a large plastic strain, the strain is gradually reduced. The plastic strain component is unchanged during the relaxation; instead, it is the elastic component of the strain that is reduced. As the stress depends entirely on the elastic component of the strain, the stress also immediately begins to decrease. Hence, we do not travel back along the “outward” path, but instead return along a shifted version of the linear elastic region: the material exhibits hysteresis. Although the initial point and the final point of Figure 10B are both characterised by zero strain, the state is quite different, as the latter state has a large plastic strain. If we continue to reduce the strain, as in Figure 10C, the stress becomes negative - a compressive stress - until, once again, we reach the yield strength limit. Note that the yield strength in tension and the yield strength in compression are not necessarily equal.

In three dimensions, an elastoplastic material is characterised by a yield surface, which is the boundary of the elastic region. The yield surface is defined in stress space (or equivalently in strain space). For isotropic materials, the yield surface is most easily defined in the space of the three principal stresses. Stresses lying inside the yield surface are associated with an elastic state. When the stress reaches the yield surface, the material yields, and further increases in strain result in the stress state moving along the yield surface. Many different shapes of yield surface are possible, and are characteristic of different kinds of materials. Yield surfaces are conveniently defined mathematically by a yield function f, such that the yield surface is the locus of points for which f = 0. A complete explanation of elastoplastic behaviour in three dimensions is not given here. We recommend that you consult a good mechanics of materials textbook.

To specify an elastoplastic material, it is necessary both to specify its elastic properties, and its plasticity. Currently Faim only allows plasticity to be defined with isotropic elastic properties. A typical specification of an elastoplastic model using n88modelgenerator will look something like this

The currently supported plastic yield criteria are described below.

The von Mises yield criterion states that yielding begins when the distortional strain-energy density attains a certain limit. It is therefore also called the distortional energy density criterion. The distortional strain-energy density is the difference between the total strain energy density and the strain energy density arising only from the part of the strain resulting in a volume change. Hence, the distortional strain-energy density is associated with that part of the strain that results in no volume change: or, in other words, the part of the strain that causes a change in shape.

The yield function for a von Mises material, in terms of the three components of the principal stress, is given by

The parameter Y is the yield strength in uniaxial tension (or compression).

To specify a von Mises material in n88modelgenerator, use the parameter plasticity. For example, --plasticity=vonmises,68 will define a von Mises material with Y = 68. The units of Y are the same as stress, so typically MPa.

The Mohr-Coulomb yield criterion is a refinement of the Tresca yield condition, which states that yielding begins when the shear stress at a point exceeds some limit. To this, the Mohr-Coulomb criterion adds a dependence on the hydrostatic stress, such that the shear stress limit is dependent on the hydrostatic stress, in particular that the shear stress limit increases with hydrostatic stress. The Mohr-Coulomb yield criterion is characterised by two values: the cohesion c and the angle of internal friction φ. For principal stresses in the order σ₁ > σ₂ > σ₃, the Mohr-Coulomb yield function is

Under uniaxial tension, the yield strength is

while under uniaxial compression, the yield strength is

Mohr-Coulomb materials are specified to n88modelgenerator using these two yield strengths. See plasticity. For example, --plasticity=mohrcoulomb,40,80 will define a Mohr-Coulomb material with Y_T = 40 and Y_C = 80. The units of Y_T and Y_C are the same as stress, so typically MPa.

The parameter material_table allows one of two standard types of material table may be chosen: a simple table, consisting of a single material definition, or a Homminga material table, which relates density to stiffness.

If you require a material table different than either of the standard ones offered, you can use a material definitions file. See below for details.

The simplest material table consists of a single material definition. Thus it corresponds to homogeneous material properties: all elements are assigned the same material properties.

Note that although only a single material definition is created, it is nevertheless still possible for the material table to have multiple entries corresponding to multiple material IDs. In this case every material ID in the table maps of course to the same single material definition. n88modelgenerator automates this: it examines the segmented input image, and generates one entry in the material table for every unique material ID (i.e. image value) present in the segmented image. Clearly, as a matter of efficiency and simplicity, if you intend to use an image for analysis with homogenous material properties, it is preferable to segment the image object to a single material ID. It is not however, necessary to do so.

Homminga (2011) introduced a model that attempts to account for varying bone strength with density. In this model, the modulii vary according to the equation

Where ρ is the CT image density. For the orthotropic case, the shear modulii follow the same scaling, and for the anisotropic case, the stress-strain matrix follows this scaling.

To use this type of material table, the data must be segmented such that the material ID of each voxel is proportional to its CT image density. As material IDs are discrete, the CT image density is therefore binned. Hence

ID_max is set with the parameter homminga_maximum_material_id, and the exponent is set with the parameter homminga_modulus_exponent.

E_max will depend on the specified material. For an isotropic material, it will be given by youngs_modulus. However orthotropic or anisotropic materials may also be scaled: the entire stress-strain matrix scales according to the above-stated law.

The method presented up to this point for defining materials in n88modelgenerator, while relatively easy, suffers from some limitations. In particular, only a single material can be specified using the options discussed so far. Furthermore, an anisotropic material can not be specified at all using the command line parameters. n88modelgenerator therefore provides an alternate method of defining materials, which is via a material definitions file. A material definitions file can be used to define any material known to Faim, any number of materials, and any material table. A material definitions file can be specified with the argument material_definitions.

The material definitions file format is straight-forward, as you can see from this example.

Here we have a material definitions file containing two materials, both of which are linear elastic materials, but with different numerical values. Besides defining the materials, we have to specify the material table, which maps materials IDs (ie. the values present in the input image) to the materials we have defined.

Here is another example material definitions file, which defines an anisotropic material.

The possible materials, and the values that must be specified for each, are listed in the following table.

n88postfaim requires that n88derivedfields has been run on the solved n88model file because it makes use of the derived fields. If any field is required for a particular table, but is not present in the n88model file, n88postfaim will simply silently skip that table.

In the usual case, running n88postfaim is as simple as,

This will generate an output file analysis.txt. The output is a report of values that are commonly of interest in finite element calculations.

The sections below will discuss a few important options that affect the output of n88postfaim. For a complete list of arguments supported by n88postfaim, see the Command Reference chapter.

Below are listed the tables provided in the output of n88postfaim.

Note that not all values or tables are generated for every model. Angular values are generated only if a rotation center is defined (see the above section), while a break-down of values per material is only performed in the case that more than one material is defined.

The first step is to open a Terminal (on Windows a Command Prompt or PowerShell). Make sure both that the Faim is on the Path, as described in How to run the solvers, and also that you have activated the Anaconda python environment for n88tools, as described in Install n88tools and vtkbone.

Assuming that you have copied the data file cube6x6x6.vti to your current working directory, the command to generate the model is

Since we haven’t specified otherwise, n88modelgenerator will select all default parameters, including default material properties (homogeneous material with isotropic Young’s modulus of 6829 MPa and Poisson’s ratio 0.3), and a displacement, applied to the top surface in the z direction, numerically equivalent to an apparent level strain of -0.01 (the minus indicates compressive strain). Refer to the section on n88modelgenerator for details on setting parameters differently if you wish to do so.

Here is the output when we run it:

We now have a file cube6x6x6_axial.n88model that can be used directly as input to the solver.

Having now defined the problem, we can solve it, and do some post-processing, all in one step with faim.

The output is:

First ensure that ParaView is set up with the plugins to read Numerics88 files, as described in Installing and using the Numerics88 plugins for ParaView. Assuming that the plugins are loaded correctly, cube6x6x6_axial.n88model can be opened directly with ParaView. When we first open this file, it looks very similar to the input data. However, if we go to the Information tab in the Pipeline Browser, we can see that there are now several Data Arrays, Displacement, ReactionForce, etc…​ as shown in Figure 17. Notice that the icon next the the array name indicates whether the data is on Points or Cells, which is equivalent in finite element terminology to Nodes and Elements, respectively.

There are many ways to visualize the data. As a simple example, choose Displacement from the drop-down menu on the toolbar that currently shows _MaterialID. Also click the Toggle Color Legend Visibility button on the toolbar (on the very left side). ParaView should now look like Figure 18. Notice that a new drop-down box appears that shows Magnitude. This occurs because the Displacements are a vector quantity. Experiment with changing this to Z, X and Y.

Although we now have the cube colored by Displacement, the object is shown in its undeformed unloaded shape. To visualize the deformations, we can use the Warp By Vector tool. There is an icon for this on the bottom Toolbar. Because we have imposed a small strain, -0.01, it will be hard to see. The Warp By Vector tool however allows us to apply a Scale Factor. We will set this to 10, so that the displacements are amplified by a factor of 10. The results are shown in Figure 19, where the compression and outward bulging of the sides are clearly seen. In this figure, we have, in addition to applying the Warp By Vector tool, re-enabled the display of the original cube6x6x6_axial.n88model, which was automatically disabled by ParaView when we connected a filter object to it. We then set its display type to Outline, resulting in the white box outline in Figure 19, which indicates the original undeformed box shape.

A typical result we want to obtain with a compression test is the object apparent stiffness. Numerically the stiffness is the net external force divided by the displacement. We can find the displacement and force on the top surface from the analysis file, cube6x6x6_axial_analysis.txt. This file is a text file, and can be opened with any text editor. Refer to the section on post-processing with n88postfaim for the complete explanation of the values available.

The displacement is obtained from table Nodal Displacements. The displacements in this table are listed according to node set.

n88modelgenerator has defined some standard node sets (for a complete list see Creating node and element sets in the chapter on Preparing Models with vtkbone.) Here the node sets that are relevant are the two surfaces to which the boundary conditions are applied, namely the top and bottom surfaces, face_z1 and face_z0, respectively. The analysis file has been generated using these two node sets. If we examine the Nodal Displacements table of cube6x6x6x_axial_analysis.txt, as shown in Figure 20, we see that the top node set has a uniform displacement in the z direction of -0.06 mm. The relevant number is circled in red. Notice also that the minimum and maximum are the same so all nodes of this set are equally displaced, as we expect, since this is in fact a boundary condition. The other node set is the nodes of the bottom surface, and it has constant displacement of 0.

Similarly to the table for displacements, table Nodal Forces shows forces, summed over the defined node sets. Here the relevant value is -2651 N, as shown in Figure 21. The apparent stiffness of the cube is therefore

Is this value reasonable? We can estimate what we should get. Completely ignoring lateral expansion and Poisson’s ratio, we can do a quick estimate that the stiffness of a homogeneous cube should be approximately,

This is close to what we obtained from the finite element model. The difference arises from the fact that the finite element model takes account of the complicated lateral expansion (bulging of the sides) that occurs in an axial test, while our simple back-of-the-envelope calculation assumes uniaxial conditions (i.e., no lateral constraint on the top and bottom boundary surfaces).

This concludes the cube tutorial. You have learned how to generate standard models from image data using n88modelgenerator, how to solve the models, how to visualize the models using ParaView, and how to do some elementary post-processing analysis.

As before, we recommend copying the data file radius_slice82.aim into a temporary working directory. The command to generate the model is

The resulting model file will be named radius_slice82_uniaxial.n88model.

Before we solve it, we are going to examine this model file using n88modelinfo. You will find that n88modelinfo is a very useful tool for getting information about the models, and for tracing back how they were created and solved. n88modelinfo can be run like this,

n88modelinfo can generate quite a bit of output. We will shortly see how to reduce this to just the information that we are interested in. For now though, let’s look at the complete output. Comments follow below the example output.

Some things to note about the output of n88modelinfo:

Now, as mentioned, that is a lot of output. Often we are interested in just part of the output. In this case, we can specify the bits we want. For example, to get just the History and the Log, we could do this,

From the previous tutorial, we know that we can use faim to solve the model and perform so post-processing on it. faim itself actually just calls three other programs to carry out the actual work. Usually it’s easiest to call faim, but in this tutorial we are going to call the individual processing programs ourselves, so we can understand what faim is doing.

The first step carried out by faim is to solve the model with n88solver_slt, which we can do with

We can see how this has modified the n88model file by running n88modelinfo again. To avoid the lengthy complete output, here we use arguments that show only the sections that have in fact changed

The output is:

Note things are worth noting here. In the Log section, we now have a record of the complete configuration of the solver (in this case, it happens to be default values). And in the Solutions section, we now have a solution (“Solution1”) with a field defined on the nodes, namely Displacements.

Faim includes a tool that can be used to evaluate the accuracy of the obtained solutions. This is discussed in the section Evaluating solution quality. It can be run like this:

The output is,

The meaning of these results is discussed in the section Evaluating solution quality; briefly we can take the number rms err/max force as a measure of the relative error. Here it is 7×10^-7, which is acceptably small.

So far our solution consists only of displacements. Many other values are of interest, for example strains and stresses. These are calculated with n88derivedfields. It can be run very simply as,

Now let’s have a look at what’s in the model file, using again n88modelinfo,

The output is:

We have a number of new solution fields (or variables). Some are defined on the nodes and some defined on the elements. For details on these fields, see n88derivedfields.

The final step performed by faim is to call n88postfaim, which does post-processing analysis. In the previous tutorial we saw that an analysis file was generated, which is a text file with a number of tables summarizing the numerical results. It is n88postfaim that generates this text file. To run it on an n88model file that contains all the solution fields is once again straightforward in most cases,

Note that unlike the previous commands, for n88postfaim we actually have to specify an output file with -o analysis.txt. (If we don’t, then the output will just be printed to the terminal.) The faim command, which combines n88solver_slt, n88derivedfields, and n88postfaim, would have automatically generated the filename radius_slice82_uniaxial_analysis.txt in this case, but that is getting a little bit lengthy, so we’ll just use analysis.txt here.

We already saw some of the values in the analysis file in the previous tutorial. For this tutorial, we are going to have a look at table Load Sharing, as shown in Figure 25. We see that this table breaks up the total load on each node set according to material ID. In this case, we have material ID 127 being cortical bone, and material ID 100 being trabecular bone. Recall that the material properties are not different; it is merely a different labelling. Here we will see the reason for doing this. Recall also that the sets defined for post-processing are the top and the bottom surfaces. (If you’re not sure of this then the command n88modelinfo --problems radius_slice82_uniaxial.n88model will show you.) Because of the alignment in the CT scanner, in this case the top surface is proximal, and the bottom surface is distal. Referring to the values circled in red in Figure 25, we conclude that on the proximal surface of this radius bone slice, the proportion of the load carried by the trabecular bone is (-806) / (-5199) = 15.5%, while on the distal surface of the slice the proportion of the load carried by the trabecular bone is (3637) / (5199) = 70.0%.

The n88model file format supports compression internally and transparently. An n88model which contains compressed data can be read with every N88 tool just like any other n88model, without the need for the user to explicitly de-compress the data. Compression is not used by default, because it does have an increased calculation overhead for reading and writing (particularly for writing). However, once we are finished with processing the model, we typically want to keep the data around for archival purposes. It is recommended to compress the n88model files for archiving. This can be done with the n88compress tool, as follows,

For this data file, enabling compression reduces its size from 438 MB to 226 MB, a reduction of 47%.

This ends the tutorial on a compression test on a radius bone slice. We were introduced to n88modelinfo, we solved a realistic model, we looked at using the residuals to evaluate the quality of the solution, we saw how to do load distribution calculations, and finally we learned how to enable compression to more efficiently store large n88model files.

We are going to apply a bending test, as described in Bending test. In order to identify the rough head of the bone as a boundary surface, we use the option bottom_surface=visible. As it implies, it will select nodes that are visible from a certain point of view. In this case, visible means “looking up from below” (effectively from an infinite distance, as parallel rays are used). Note that the distal head of the bone is in the “bottom” of the image, that is at smaller values of z. With this option, you can imagine that we might “see” parts of the object that are far away if they are not obscured by the foreground part of the image. We don’t want these distant parts of the object to form part of our boundary set, so when we use the bottom_surface=visible option, we nearly always want to set a depth limit as well using bottom_surface_maximum_depth. In this case a suitable limit is 15 (units are the same as the image units, which are millimeters in this case). In general there may be a bit of experimentation to get the best depth limit.

Since we have several parameters that we want to specify, it will be more convenient to create a configuration file for n88modelgenerator than to specify them individually on the command line, as we have done previously. To do this, open a text editor, and create the following file. Save it to your working directory as radius_bending_test.conf.

Lines starting with # are comment lines: any line starting with # is ignored.

We can now run n88modelgenerator by specifying only the configuration file to use:

Because this is a non-trivial boundary set, we really want to visualize it to make sure it is reasonable and that it corresponds to what we expect. The best way to do this is with the tool n88extractsets. This is the command to use:

Here “--constraints” means both boundary conditions and applied loads. (We can leave off this option, but then we will get even more output files which in the present case we are not interested in.)

The output looks like this:

We see that four files are generated, one for each constraint. These are VTK “PolyData” files, which consist of a set of “Vertices” (points), one for each node in the boundary condition set. n88modelgenerator has separated the “fixed” constraints, which are degrees of freedom that are fixed at zero, and the “displacement” constraints, which are degrees of freedom that are fixed at some non-zero value. This distinction is sometimes useful, although the solver does not require it. In this case the nodes belonging to the constraint or boundary condition have x and y fixed at zero, and z at non-zero values. Thus the “fixed” and “displacement” sets will consist of the same set of nodes, but with different senses and values. We are going to open the file radius82_bending_constraint_bottom_displacement.vtp displacement and render it together with the contoured input image. This is shown in Figure 27. Each individual node belonging to the boundary condition is shown as a blue dot.

This looks good so we can go ahead and solve this model with faim. Expect this to take a while - go grab a coffee.

As this is a very large file, once solved we should enable compression as in the previous tutorial. Compression can take a rather long time on a file this large. Remember that it is never necessary to uncompress an n88model file before reading it.

A comparison of the original input image with the solved model is shown in Figure 28. Here we have used a Warp By Vector filter with a Scale Factor of 15 to make the distortion more evident to the eye.

Another useful way to render this data is coloured by von Mises stress. See VonMisesStress. This is shown in Figure 29.

Showing the von Mises stress as in Figure 29 is informative, but only indicates the stress on the visible surfaces. Often we will want to look at the stresses within the volume. One way to do this is with a Slice filter. This is shown in Figure 30. Without going into a great deal of detail, you can see the filters used in this rendering in the Pipeline Browser.

The interesting quantity in the analysis file radius82_bending_analysis.txt (see the section on post-processing with n88postfaim) is the torque around the bending axis, which in this case is in the y direction. Refer to Figure 8. This quantity is given as the total of Ty in the following table from radius82_bending_analysis.txt. By default, this is calculated around the center of mass of the object, although you can re-run n88postfaim and specify an argument for the argument rotation_center if you want the torques calculated around an axis passing through a different point. Because the head of the bone is “down” in the sense of smaller z, face_z0 is the head of the bone, and face_z1 is the cut face of the shaft. In any case, we can see that the torques balance. The numerical value of the y compenent of the torque is 1997 N·mm (see A note about units), or 1.997 N·m.

Because nonlinear models take substantially longer to solve than linear models, for this tutorial we are going to use the reduced-resolution image file radius164.aim. However, if you have the processing power and the time to obtain solutions, and you want to use the original resolution file radius82.aim, you can of course do so.

Creating a model with elasto-plastic material properties with n88modelgenerator is straight-forward: one simply has to set the plasticity option. For example, the following configuration file is similar to the one in the last tutorial, but we have added plasticity = vonmises,10. This means that we want a von Mises yield criterion (see von Mises yield criterion) with a yield strength of 10 MPa. There is no scientific justification for choosing this numerical value for the yield strength: it was chosen simply because it produces interesting results for this tutorial.

As before, we generate the model with

We can verify that the resulting model file has elasto-plastic material definitions using the handy n88modelinfo command,

The output is

We can solve the model with faim,

Let us look into the file to see what fields are present for a solved elasto-plastic model. The command to list the solutions is, as before,

this gives an output of

Everything here is as before for linear solutions, with the addition of a new field PlasticStrain. PlasticStrain is the irreversible strain that arises from exceeding the elastic limits as defined by the yield surface. See PlasticStrain in the section Calculating additional solution fields with n88derivedfields. You may also note that variables defined on gauss points is a new category. In essence, this is a higher accuracy version of variables defined on elements, and it is required for accurately processing elastoplastic models.

We can use ParaView to open the solved model file radius_bending_elastoplastic.n88model. In Figure 33 we’ve rendered the surface of the model such that unyielded elements are shown in grey, while yielded elements are shown in red. For elasto-plastic materials, there is a clear distinction between elements that remain within the elastic limits and therefore have zero plastic strain, and those which have exceeded the elastic limits and have non-zero plastic strain. (where non-zero means that at least one of the 6 components xx, yy, zz, yz, zx and xy of the plastic strain is non-zero). In order to generate the rendering of Figure 33 we set the Coloring to the Magnitude of PlasticStrain. You can find these in drop-down boxes in the toolbar in ParaView, as shown in Figure 31. Next we will manually set a tiny range for the color scale, thus ensuring that practically every element that is yielded is shown in fully saturated red. If we set a symmetric scale around zero, then unyielded elements, which have a PlasticStrain magnitude of zero, will be drawn in a neutral gray. To set the color scale, click Rescale to Custom Data Range, which is just to the left of the drop-down box for PlasticStrain, as shown in Figure 32. For this rendering I used a range of -0.00001 to 0.00001 . Note that in Paraview, the Magnitude of the PlasticStrain is simply the square root of the sum of the squares of the 6 components of plastic strain. This quantity has no physical interpretation, but it is useful for distinguishing yielded and unyielded elements, as here.

What we observe in Figure 33 rendering is that there are essentially three regions of plastic deformation in this bone subject to a bending force. They are,

You can’t actually tell from Figure 33 what is yielding in compression and what in tension, but with a bit more manipulation, we can generate an image that does show this. Since we know a priori that the tension and compression on the sides of the bone shaft will be predominantly in the z direction, we can change the component of YieldStrain from Magnitude to ZZ. (You may need to manually re-adjust the color scale, as ParaView will automatically rescale to the data range when you change the variable for the color.) Elements yielded in compression now show as blue, as in Figure 34. We have in addition used a Slice filter, as in the previous tutorial. This allows us to see that the yielded region is rather thick.

In comparison with the linear case, the y component of the torque has dropped about 1.5%, from 1.997 N·m to 1.968 N·m; hence a decrease in stiffness of 1.5%. In conclusion, this amount of plasticity has not significantly affected the results.

There are always more subtle considerations for nonlinear calculations than for linear calculations. One aspect of elasto-plastic calculations is that they exhibit hysteresis, or irreversible changes to strain. In this way, the solutions depend in fact, not only on the current applied state, but in principle on the whole time-dependent history.

Whether this is significant or not in practice is something that can be investigated experimentally - or rather computationally in our case. In the present case, we already know that there is only a 1.5% difference in numerical results between the linear model and the elastoplastic model. So this more careful method of applying the load in the nonlinear case is very unlikely to result in significant differences in the result. This tutorial is therefore more of a demonstration of principle, and an introduction to the possibility of modifying already-solved n88model files.

To apply the bending in increments, we will take advantage of the fact that n88modelgenerator can accept an existing n88model file as input. In this case, it will preserve the grid, i.e. the points and elements of the model, but generate new boundary conditions and/or material definitions. Crucially, it also preserves an existing solution, which is then automatically used as the starting value when the solver is subsequently called.

We are first going to make a configuration file containing all the parameters of n88modegenerator that are the same for each increment of bending.

Now we are going to call n88modegenerator and faim multiple times, in steps of 0.1º, from bending angle 0.6º up to bending angle 1.0º. 0.6º is chosen as a the starting angle because it happens to be just after the onset of nonlinear behaviour. In general, it takes a bit of experimentation to find the minimum applied test condition for which some elements start to exceed the linear region.

Here are the commands. You may want to put these into a script file.

Once we have run all of these commands, we can verify that the final data file, incremental_bending_1.0.n88model, does indeed show the history of being processed multiple times:

This gives an output of:

As usual, the details can be recovered with the --log options of n88modelinfo, but the output is lengthy, so we will not show it here.

To perform the post-processing on the final result, we call,

Numerically, the results have changed little in this case, as expected. We now have a y-component of torque of 1.972 N·m, which compares with 1.997 N·m in the linear case, and 1.968 N·m for the “all at once” elastoplastic case.

This tutorial is the final tutorial in the series presenting the Faim standard workflow, with models generated by n88modelgenerator. Those interested in generating custom models will want to continue on to the tutorials presenting vtkbone. If you intend for the moment to use only n88modelgenerator for generating your models, then you can stop here.

It is not necessary to know Python in order to write and modify scripts using vtkbone, as very often you can make obvious modifications to existing scripts. However, having a basic familiarity can make you more comfortable with the following scripting. As Python is very widely used, there are any number of excellent introductory books. A suggested place to start is the on-line tutorial at http://docs.python.org/tutorial/ .

One feature of Python which we will use that is not a part of the standard Python distribution (but is distributed with Numerics88 software) is Numpy, which efficiently handles numerical arrays in Python. Again, prior knowledge of Numpy is not required, but will help. Documentation for Numpy can be found at http://docs.scipy.org/doc/ . We recommend https://docs.scipy.org/doc/numpy-dev/user/quickstart.html as an introduction.

Finally, there is VTK . VTK is quite complex and rather more difficult to learn than either Python or Numpy, but fortunately it is not necessary to know very much VTK. We will be using a mixture of custom vtkbone objects and standard VTK objects. When one of these comes up in the following tutorials, you may want to have a glance at the corresponding API documentation for that object. This is explained in vtkbone API documentation.

First create a working directory. Using your favourite text editor, create a file apply_compression.py.

Start by adding these lines, which import all the Python modules that we will be using:

Once you save this, you already have a script that you can run! In fact, you should run it, with just these lines in it - even though it does nothing yet, because it will verify that Python can successfully locate VTK and vtkbone. The script can be run like this (make sure you’ve saved apply_compression.py),

You should get as output . . . nothing whatsoever. If you get an error like “ImportError: No module named vtk”, then it is likely that the n88tools is not installed, or that the anaconda python environment is not activated. Refer to Install n88tools and vtkbone.

We recommend also having a version number in your script. Increment it whenever you change the script (not every time you add a line, but if you modify the finished script later so that the resulting model is somewhat different, then increment the version number.) We’ll see soon how we can ensure that the script version number gets written into the log of the n88model file. Although we’re not accountants, we highly recommend having an audit trail that will allow you to know exactly how a given model was generated. (Obviously, you should archive copies of every version of your script, not just the most recent one. You should also floss daily, which we’re sure that you do.) For now, just add a version number to the script with this line. We are on version 3, since this tutorial has been in the manual for quite some time, and it has been tweaked a couple of times.

Now, we’re going to want to specify the input image file on the command line, because this will allow us to use the script on any input file, not just one hard-coded by name into the script. We’re also going to want to inform the user that an input file name is required if they don’t supply one. (“The user” is most likely you yourself, but don’t assume that next week you’ll remember how you intended to use the script you write this week.) The python module argparse is really handy for this. The following lines will do the desired thing:

So now when we run it as before, this time we do get some output:

Progress!

Let’s see what happens if we run it with the -h flag. Note that argparse has automatically added a -h flag, the purpose of which is to get help about the command.

The input file is cube6x6x6.vti and as usual you can find it in the data folder of your distribution. We’ve seen it before in the first tutorial. Copy it to your working directory.

As the data is in VTK’s .vti format (properly referred to as “VTK XML Image Data file format”), we’ll need a vtkXMLImageDataReader object to read it. That will look like this:

vtkXMLImageDataReader is a type of VTK “filter” (specifically a “source”). What we’ve done is very typical for VTK filters and we’ll see this sequence of steps for employing a VTK filter over and over again. It is,

In fact, this isn’t the only way to use VTK filter objects, but it is the simplest and it is the procedure we will nearly always employ in our scripts.

Let’s provide some feedback at this point so we have some indication that the data was correctly read. While we’re at it, let’s do some error checking to stop things right here if that data wasn’t able to be read,

Now running it should result in,

We could assume that the input image is properly segmented and just proceed with generating a mesh. Since this is an introductory tutorial perhaps we should. However, it happens not infrequently for whatever reason that a finite element model is attempted on an image which does not consist of a single connected object. This can cause all sorts of trouble. Consider this: if there are disconnected floating bits in the model, where is the proper place for the solver to put them in the solution? Sometimes the solution of such models will be very slow, other times the model will solve OK, but some of the post-processing numbers may appear strange if you don’t realize that they are calculated over the entire model, including bits which are subject to no constraints.

So to avoid potential grief, this is a good place to drop in a vtkboneImageConnectivityFilter . It will strip away (actually zero-out) any parts of the input image that aren’t in some way connected to the largest object in the image. It is a safe and prudent step in the generation of a well-defined finite element model.

The filter object for generating an FE mesh from an input image is vtkboneImageToMesh . The FE mesh represents the geometry of the model, and consists of elements and nodes (concretely, it will be an object of type vtkUnstructuredGrid ). Note that in VTK terminology, elements will be called Cells, and nodes will be called Points. Most of the time we don’t need to set any parameters to vtkboneImageToMesh except the input.

Now our output is beginning to get interesting,

Quick check: 6³ = 216 and 7³ = 343, so we’re good.

Now we need to define some material properties. We can do this by creating an instance of some specific derived class of vtkboneMaterial . The one of these is vtkboneLinearIsotropicMaterial . Let’s make one of them,

The other thing we must do is to make a material table, which associates the values in the input image (segmentation IDs or material IDs) with defined materials. Our image file cube6x6x6.vti has the uniform value of 1 on every voxel, so we need to make a material table that associates 1 with our material defined above. But perhaps we might want to run this script someday on an input file that uses different values (127 is very common for historical reasons). We could write here some fancy scripting that would examine the input data and produce a list of all the different values in the image, then map all those values to our material (since we’re only using one material here). This sounds a bit complicated, but fortunately, there is a vtkbone object, vtkboneGenerateHomogeneousMaterialTable , that does exactly this for us:

To apply boundary conditions appropriate to a compression test, we can make use of the object vtkboneApplyCompressionTest , which will do practically all the work for us. vtkboneApplyCompressionTest requires as inputs both the mesh and the material table that we defined above.

The output of vtkboneApplyCompressionTest is a vtkboneFiniteElementModel object, which is a complete representation of the problem we want to solve with finite element analysis. It remains to write out this model in n88model file format.

That’s it. If everything works correctly, the output should look like this:

We promised at the beginning to record the script version number in the file, and so far we haven’t done this. First we’ll have a look at what is currently in the history field.

This is useful information so instead of replacing it, we can just add to it. To do this, add the following line to the script, just before writing out the file:

Now run the script again, and this time when we run n88modelinfo, we get:

Good, now we’ll always know that the cube6x6x6_custom.n88model file was created with version 2 of the script (which called vtkboneApplyCompressionTest from vtkbone version 6.0).

We can also have a look at the log, which gives more detailed information:

There is already quite a lot of information here, provided automatically by vtkboneApplyCompressionTest . In particular we have all the parameters of vtkboneApplyCompressionTest , which were all set to default values because we didn’t change any of them. We could, if we wanted, add even more information to the log, using the method model.AppendLog, just as we used model.AppendHistory. In this case, there is nothing lacking that we can’t query. For example, although no material information appears in the log, we can determine it with n88modelinfo,

You can now proceed if you want to solve and analyze the model exactly as in An introductory tutorial: compressing a solid cube.

This completes the tutorial on building the 6×6×6 cube compression model using vtkbone. We’ve seen that it takes considerably more work to write a script than to call n88modelgenerator, so whenever it suffices, it is preferable to use n88modelgenerator. However, although in this case we’ve done nothing that n88modelgenerator couldn’t have done, we hope that you can imagine that now that we have a script, that we could potentially modify it in any number of ways beyond what would be possible with n88modelgenerator. In subsequent tutorials, we’ll explore these possibilities.

As in the previous tutorial, we start the python script by importing the python modules that we will need. This time we will also need the numpy module, as we are going to do some array manipulation.

As before, we’ll also set a script version number:

As before, we use the argparse module, which makes it quite simple to define a couple of arguments for running the script, one of which, the applied force, will have a numerical value.

The beam is going to have a very simple geometry, one which we generate algorithmically without too much trouble. This is also going to serve as a brief introduction to numpy arrays. Because the analytic solution for the cantilever is essentially a 2D problem, we are going to use just a beam width of just a single element. To scale up to a realistic width, if desired, one can just scale all the forces. We are going to create a beam of 200×1×20 elements, with a cubic element of side length 0.5mm. Thus the beam has dimensions of length 100mm, height 10mm, and width 0.5mm (but again, everything scales linearly with the width, so the obtained solution can be used for arbitrary widths).

Here is how we create a 3-dimensional array data of size 200×1×20:

We need to convert this from a numpy array, which is natural to Python, to a VTK data array:

There are a couple of subtle points here. The first is that VTK is actually going to want the flattened data, meaning as one long 1-dimensional array, rather than as a 3-dimensional array. (It will keep track of the dimensionality separately, as we will see below.) Hence the use of ravel function, which is a numpy function that returns the input array as a one-dimensional array (without copying the underlying data if possible). The other aspect to note is that we specify deep=1, which specifies that the VTK array should create its own copy of the underlying data, instead of just referring to the numpy data. This is less memory efficient, but there are potential pitfalls with sharing the data, so we’ll go with the more robust method of making a copy.

Now that we have the data, we’re going to create a vtkImageData object which adds information such as the dimensions and the spacing.

We’ve set the spacing to 0.5 . In this tutorial, we are going to use units of millimetres. Units are not specified in the n88model file. Instead we will just have to ensure that all our units are consistent. See A note about units.

We can at this point write out the vtkImageData as a .vti file, and open it with ParaView to ensure that we have obtained the desired result. It should ressemble Figure 35. As a beam, it looks absurdly thin, but to emphasise once again, we are solving essentially a 2D problem; the numerical results obtained can be scaled to any beam width.

Now that we have an image, as in the previous tutorial, we convert it to a geometric mesh:

We also define a material (steel) and a material table.

In contrast to the previous tutorial, we won’t be using one of the standard test generation filters. Rather we will use the base test generation filter vtkboneApplyTestBase , which creates a vtkboneFiniteElementModel object without any boundary conditions.

Before we add displacement boundary conditions and applied loads, let’s save the model so we can examine the results as we build up our script. Of course the model isn’t complete yet, but that doesn’t prevent us from saving and examining it.

To run the script will require both a numerical value for the applied force and an output file name. For now use these values:

We are choosing to use a negative applied force, because it is more natural to think of a downward force applied to a cantilever.

We want to fix one end of our beam, to weld it to the (virtual) wall. Usually the step requiring the most work in defining boundary condition is identifying the set of nodes that constitute the boundary surface we are interested in. In this case, vtkboneApplyTestBase has done the work for us, because although it doesn’t define any boundary conditions, it does define node (and element) sets on all the faces. (To be more precise, it defines node and element sets consisting of any nodes/elements located at the x,y and z extents of the model.) This is discussed in Creating node and element sets. The names of these sets are listed there, but if you’ve run the script up to this point, you can also use n88modelinfo to query the set names:

There are also corresponding element sets, which can be listed with the --element_sets option to n88modelinfo.

Node set face_x0 is the one to which we are interested in applying a fixed boundary condition. If you refer to the the API documentation for vtkboneFiniteElementModel, you will see that there are a number of methods for defining boundary conditions. One possibility is to use the method FixNodes, which will fix the specified node set in all senses (directions). Only one line of code would be required:

However, as we are trying to obtain a result as close as possible to the analytic 2D solution, it is actually better to not limit the freedom of motion in the y direction on the boundary condition. This way, no y direction forces will arise. Were there to be any forces in the y direction, this would cause some essentially 3D effect on the x-z plane by the mechanism of Poisson’s ratio. So instead of using FixNodes, we will use the following slightly more verbose code

If you now run the script, you can then use n88modelinfo to verify that the boundary condition has indeed been generated:

At the free end of the beam, we are going to apply a load. Creating an applied load is very similar to creating a displacement boundary condition. vtkboneFiniteElementModel provides a number of methods for this purpose. There are some important differences,

Here is the code to create an applied load constraint, named end_force to the pre-defined element set face_x1. corresponding to the end face of the beam.

As discussed in Convergence measure there are different possibilites for the convergence measure, which is used by the solver to determine when the solution has iterated enough to obtain an acceptable solution. Generally, the best convergence measure to use for linear problems is convergence set. This convergence measure can only be used if we define a convergence set within the model file. A convergence set is something that we choose, and we generally want to choose the quantity that we are most interested in calculating. For example, for a cantilever, this quantity could be the displacement of the tip of the cantilever. You might imagine that this would be difficult to specify. But in most cases, this quantity of interest is the complement of some boundary condition or applied force. In this case, the displacement is the complement of the force applied to the tip. There is a method to generate a convergence set from a boundary condition or from an applied force, and it is very simple to use:

This method notices that the constraint named end_force, which we created previously, is an applied force, and that it is applied to the faces of certain elements. It then constructs the compliment, which is the displacement, averaged over all nodes located on those element faces. It is this average displacement of the tip that the solver will watch to determine when its value ceases to change meaningfully. At that point the solver is finished.

There is one more step required, and that is to identify the sets to be used for post-processing calculations by n88postfaim. We are interested in values on the sets face_x0 and face_x1. The magic lines of code are:

After running the script, we can verify that these post-processing sets are in fact specified in the n88model file:

As discussed previously, it is good practice to document the model creation in both the file History and the file Log. Here we add these lines (again they must come before the code that saves the model):

If we run the script to generate the n88model file, we can at any later time get a brief description of the model and determine what versions were used (including what version of vtkbone).

The model is quickly solved with faim:

If we open the solved file cantilever.n88model in ParaView, we can observe the beam deflection as shown in Figure 36. Here the deflection is amplified by 10 so as to be more evident. The figure is colored by the xx component of strain, with red being tensile and blue being compressive. This is easily done in ParaView, simply by selecting the relevant quantities from the drop-down boxes in the toolbar, as in Figure 37.

Now let’s compare quantitatively with the analytical result for deflection of a cantilever beam as found in any Mechanics of Solids textbook or on wikipedia at https://en.wikipedia.org/wiki/Deflection_%28engineering%29. We have that the deflection δ of the end of a cantilever beam is

where

For a beam with a rectangular cross-section \(b \times h\), the area moment of inertia is

so that we have

You can observe this deflection in ParaView, and we can get a precise number by looking in the beam_analysis.txt file, in the table Nodal Displacements. In this case the result, 0.4018, is very close to 0.4 .

This concludes the cantilever beam tutorial. We learned how to create custom models using vtkbone that are beyond the standard types of models that n88modelgenerator can generate. We were also introduced to a new type of constraint, which was the applied load.

The cantilever problem has a known analytic solution even in the case of nonlinear elastoplastic material properties. This solution is due to Timoshenko and Gere (1997).

Before building and solving some models, we will consider the analytic solution, because it will indicate what range of test loads are both interesting and reasonable, which is a non-trivial question for nonlinear problems.

The applied load at which the cantilever begins to yield, F_y, is

where σ_y is the material yield strength (i.e. the value of uniaxial stress at which the material starts to yield), and all other symbols are as in the previous section.

For our current model, for which we will use σ_y = 120 MPa, this evaluates to

You will notice that this value corresponds (not entirely coincidently) with the applied force used for the linear case in the previous section, for which we had calculated a deflection of 0.4 mm.

In the yielding region, the beam deflection is

This expression is only valid for 1 ≤ F / F_y ≤ 3/2; the cantilever will not withstand any force larger than this. For the current model this range corresponds to between 10 N and 15 N.

The model generation script is named create_cantilever_ep.py. It differs from the script in the previous tutorial, create_cantilever.py, only in the definition of the material. This becomes

For a von Mises elastoplastic material, we have exactly one new parameter as compared with a linear material, and that is the yield strength, here set to 120 MPa. The exact parameters required depend on the material. For example a Mohr Coulomb material would have two additional parameters, which are set with the method SetYieldStrengths(Y_tension,Y_compression). As always, for exact usage refer to the vtkbone API documentation ( vtkboneVonMisesIsotropicMaterial and vtkboneMohrCoulombIsotropicMaterial ).

As in the linear case, we generate the model with

And we can quickly verify that this model uses a von Mises material with

The output is

The cantilever is, for various reasons, a particularly difficult model for a finite element solver. When combined with nonlinear material properties, in order to get the most accurate results, it is necessary to tighten the convergence tolerance. Thus for the following results, we ran the solver as follows,

A solution with applied force of 14.5 N is shown in Figure 39. In comparison with the linear case of Figure 36, we see that the strain is much more concentrated near the fixed end of the beam.

By generating and solving models for various loads in the plastic range, we can generate table Elastoplastic cantilever results. To do this, first we create a number of sequence of models

You probably want to put all these commands in a script or batch file.

And then we solve them all

As we can see from the number of elements exceeding the plastic limit (the last column of the table), FAIM does indeed show plastic onset at 10 N load. The data from table Elastoplastic cantilever results are visualized in Figure 40. Agreement is very good until the onset of beam failure approaches. Then there begins to be a some divergence between FAIM and the analytic results. Better agreement can be obtained with a smaller element size, but fundamentally near the singularity of failure, there are very rapid changes in plastic strain through the cross-section, and these cannot be well modelled by numerical discretisation. However, away from this singularity, numerical results are good.

Up to this point, we have solved a sequence of nonlinear models at different loads without giving a lot of thought about how these loads are applied. In fact we have solved each load point "from scratch". From a theoretical point of view, it would be better to start from the greatest load for which the model stays in the linear range, and then incrementally add small amounts of load, obtaining a sequence of slowly changing nonlinear solutions, each of which is used as the starting point for the solver for the subsequent load. This is discussed in Obtaining accurate nonlinear solutions by progressively applying loads. The question is, does this actually matter? Here we will investigate this and discover that the answer is "it depends".

The first thing we need is a script to take an existing solved n88model file, and modify the applied force. You can find this script, named modify_force.py, in the cantilever subdirectory of the examples. Since you are getting by now pretty good at scripts, here are all the preliminaries, up to reading in an existing n88 model file:

Now, we could modify the numerical values of the force constraint, but from a programming view, it is easier to just delete it and re-create it.

That’s the meat of it. Now we just finish up with the usual entires to the history and the log, and write out the modified file.

Now to run this, we first create and solve a nonlinear model as we did above, with an applied force that we know, from the previous analysis, is just at the limit of linear behaviour.

As previously, we can extract from the analysis file just the numerical value for the deflection, like this

which returns -4.020E-01.

Now we use our modify_force.py script to increase to applied force to -10.5 .

And we solve it and extract the new deflection.

What we want to see, is that faim starts with the following two status lines:

And further on, once the solver solver_sp itself has been launched, we want to verify that we see the following status output:

Then, rinse and repeat until the desired final load is reached:

and so on . . .

We can thus add another column to Elastoplastic cantilever results, and compare the “all-at-once” nonlinear solutions with the incremental solutions. Elastoplastic cantilever incremental results shows just these two columns. The results are practically indistinguishable. So can we therefore conclude that carefully applying an incremental load is pointless extra work? Not quite. We are not at the end of the story yet. Let us take our solution at a load of 14.5, and from there decrement the applied load back down to the linear region. Programmatically, this looks much like before:

and so on . . .

We can tabulate all the deflections at once,

These additional load decreasing points are tabulated in Elastoplastic cantilever results for decreasing loads. These are plotted together with the increasing load points up to 14.5 in Figure 41. Now we see that that there is in fact a difference. Indeed we recognize the classic hysteresis curve. (Note that compared with Figure 40, we have swapped the axes, just so that it looks more conventional.)

This concludes the tutorial on a cantilever with elastoplastic material properties. We learned how to use elastoplastic materials in vtkbone scripts, and also how to apply loads in incremental steps to nonlinear models by modifying an existing n88model file. By carefully applying incremental loading to the cantilever, we determined that results were essentailly unaffected as compared with an “all at once” solution for the initial deflection of the beam, but that a residual deflection remains during the subsequent unloading of the beam. This residual deflection, which is the hysteresis, can only be observed by applying the load incrementally.

This example is a finite element model of an important mechanical test for orthopaedic biomechanics. This tutorial is based on the work of Stadelmann et al (2012). Physical tests were done in which a volume of bone substitute was constrained with a ring and the screw was pulled out. Here we construct a finite element model corresponding to this physical test.

The geometry of the model is in the file sawbone.aim in the data directory of the distribution. It consists of a volume of bone substitute into which a screw has been embedded. This is shown in Figure 42. The data have been segmented such that bone substitute is labelled with ID 127 (shown in beige), and the screw with ID 90 (shown in blue). The mechanical test is to find the force required to pull out the screw.

Using your favourite text editor, create a file screwpullout.py.

By now you should be pretty familiar with importing the necessary Python modules:

As before, we’ll also set a script version number:

Since generating a model from a large micro-CT scan can take a lot of time, it is convenient to print status messages as the script proceeds, and to have these messages time stamped, so that we can see what steps potentially should be optimized (or where the script fails). The following function will take one or more strings and print them out one per line, with a time stamp on the first, and indenting the subsequent lines.

This model is going to be quite complicated. To keep things orderly, we will put all our parameters in a configuration file. The configuration file is a record of exactly which parameters were used to generate a particular model. We’re going to write the configuration file itself in Python code, so that we don’t have to write a file parser to read it. Here is what the configuration file (screwpullout.conf) looks like:

Now we add to our script the following:

All that looks like this:

Now when we run our script, we also have to specify the name of the configuration file:

We will use a vtkboneAIMReader to read the data file, which is in the Scanco AIM file format. This file contains a segmented image. As well, we want to add some additional lines to report our progress, and print some values as a sanity check.

The output so far, when we run the Python script, will look something like this:

As in the tutorial Compressing a cube revisited using vtkbone, we need to convert the input CT image data to a vtkUnstructuredGrid object representing the geometry. Again, we will use a vtkboneImageToMesh filter for this purpose.

Once again, we repeat the steps of the tutorial Compressing a cube revisited using vtkbone, and we generate a material table and two materials: bone and screw.

In the tutorial, Deflection of a cantilevered beam; adding custom boundary conditions and loads, we used vtkboneApplyTestBase to combine the geometry with the model table into a vtkboneFiniteElementModel object. vtkboneApplyTestBase didn’t create any displacement boundary conditions or applied loads for us, but it did define some standard sets that we subsequently used to define our displacement boundary conditions and applied load. This time, those standard sets won’t we very useful for us. So instead of vtkboneApplyTestBase , we’ll use an even more elementary model-creating filter called vtkboneFiniteElementModelGenerator , which creates a vtkboneFiniteElementModel object without any pre-defined sets.

Now that we have a vtkboneFiniteElementModel object, we can save it. Of course, the model isn’t finished yet, but just as in the previous tutorial, we can still save it and examine it as we go.

By now you’ll be familiar with adding History and Log fields to help us later trace back the details of the creation of this model. Let’s add all this information (before the code to save the file!):

Now run the script. As a reminder, this is done with

As shown in previous tutorials, we can retrieve the History and Log fields from the n88model file at any time with n88modelinfo:

In order to properly apply the desired constraints, we will need to know where the screw is. We can do this algorithmically, by calculating the principal axes of inertia of the screw using vtkboneTensorOfInertia . We are going to assume that the principal axis closest to the z axis is the screw axis. This won’t be exactly true. Consider for example if the screw is cut off unevenly, then the principal axis of inertia won’t be exactly the same as the screw axis, but for our purposes, it is good enough.

Now we run it, and this part of the script generates output like this:

Comparing by eye with the rendering in ParaView, we can determine that these values are about right, so we give ourselves a gold star in advanced math and carry on.

We want to add a boundary condition to the screw where it intersects the boundary of the model. All these nodes lie on a plane, so they are easily found with the function AddNodesOnPlane provided in vtkboneNodeSetsByGeometry . Note that these nodes get associated with the FE model object, and we have to assign a name, in this case “screw_top”, so that we can refer later to this set.

We want to be able to inspect this result, and ensure that we have in fact selected the nodes on the top of the screw. To do this, run the script to update the model file. Now, as in the tutorial Bending test of a radius bone with an uneven surface, we are going to use the tool n88extractsets to extract and visualize the node set.

The file sawbone_node_set_screw_top.vtp can be opened an visualized with ParaView. The nodes are shown in Figure 43, exactly on the end of the screw. This is what we want.

Having verified that we have found the set of nodes that we want, we can use the method ApplyBoundaryCondition of the vtkboneFiniteElementModel object to create a boundary condition. Note that we also give this the name “screw_top_displacement”.

This boundary condition is much more difficult, because we want to find nodes lying approximately in a ring on a rough surface. The approach we will take is to successively whittle down candidate nodes by applying a sequence of criteria, until we have finally just the nodes we want. These are the criteria:

VTK uses the concept of implicit functions which can be used to select points and cells within or without a certain geometry. Here we are going to use two vtkCylinder objects, and combine them with a difference operation with a vtkImplicitBoolean . This will result in a hollow cylinder. VTK has a filter object, vtkExtractGeometry , that will select points and cells based on a specified implicit function.

One complication is that set of Points in the output of vtkExtractGeometry are different than the original Points, and so the Point IDs will be different. However, we need to know the original Point IDs in order to be able to apply boundary conditions to the FE model. The solution to this is to use a concept in VTK called “Pedigree ID”. These can be accessed with the method GetPedigreeIds from the PointData.

Now, this is a good place to see what we have so far for a selection. There are two possibilites: we could write out geometry_in_ring with a vtkXMLUnstructuredGridWriter filter, or we could add the element IDs as a set to the model. Let’s do the latter:

Again, run the script to update the n88model file. We can now extract this element set from it as follows:

Figure 44 shows this selection (i.e. the file sawbone_element_set_in_ring.vtp), together with the screw elements, shown for context. (We used a Threshold filter in ParaView to show just the screw cells from sawbone.n88model.) Figure 45 is another view of the same data. This view is more suited to inspecting the ring geometry of our selection. To construct this view, first hit the +z button on the toolbar in order to view exactly along the z axis. Then from View Settings under the Edit menu, select Use Parallel Projection. Finally, change the Representation of both data sets to Surface With Edges in order to show individual Cells.

To select elements based on depth, as measured from the highest node in our set, we apply a similar approach as we did to select elements within the correct radii. We use a VTK implicit function, in this case a vtkBox , to specify the geometrical volume to which we want to limit our element set.

Now just as before, we can create a corresponding element set for the purposes of inspection:

Run the script, and use n88extractsets to extract the element sets, and it will create a new file, sawbone_element_set_depth_filtered_ring.vtp.

Figure 46 shows the remaining elements after filtering by depth.

We can use the function FindNodesOnVisibleSurface from vtkboneNodeSetsByGeometry for this purpose. We need to pass it a normal vector for the viewing direction, as well as an empty vtkIdTypeArray for storing the result.

Again, run the script to update the n88model file. This time we want to extract the node sets

The resulting ring-shaped cloud of points (file sawbone_node_set_ring_top_visible.vtp) is shown in Figure 47. These are the nodes to which we will apply our fixed boundary condition.

Now that we have a named node set for this constraint, it is easy to create the desired boundary condition, which we will name “fixed_bone_surface”:

As we did in Tutorial: Deflection of a cantilevered beam; adding custom boundary conditions and loads, we will add a convergence set. See Convergence measure for discussion.

This method notices that the constraint named screw_top_displacement, which we created previously, is a boundary condition. It then constructs the compliment, which is the force summed over all nodes subject to the boundary conditions. It is this force that the solver will watch to determine when its value ceases to change meaningfully. At that point the solver is finished.

Just as we saw in the tutorial Deflection of a cantilevered beam; adding custom boundary conditions and loads, there remains one more step to complete the n88model file, and that is to identify the sets to be used for post-processing calculations by n88postfaim. We are interested in values calculated on the sets “ring_top_visible” and “screw_top”. However, when we have generated our own custom sets, as we have here, we must be careful about providing both node sets and the corresponding element sets to n88postfaim. If you review how we created the sets for generating boundary conditions, you will observe that we only generated the node sets. Fortunately, there is a method GetAssociatedElementsFromNodeSet of vtkboneFiniteElementModel that will identify the corresponding element set given a node set.

After running the script, we can verify that these post-processing sets are in fact specified in the n88model file,

We now have a complete n88model file. If you want, you can proceed to solve it using faim, as in previous tutorials. We leave it as an exercise to identify, from the analysis file, the net forces between the screw top and the constrained ring of bone.

This concludes the screw pull-out tutorial. In this tutorial we learned how to generate node sets based on complex criteria, to which we can then apply boundary conditions. This is the final tutorial in the series on learning vtkbone.

We will write the script in the following sequence:

Start by opening a new script file titled “_modelgenerator_registered_fe.py_” (or another name that describes the purpose of the script).

Now we will import the packages that are needed for this script. Any additional packages that you wish to use may be added below this list:

Let’s give the script a quick test run to make sure everything is installed correctly:

Running this script should result in no output, so any output will indicate an error with the import commands. Most often, any errors that occur at this point are likely triggered by a package that is not installed in your Python environment, or your Python environment is not activated.

A useful addition to any script, particularly a script with many sequential steps, is a utility function that can output messages (i.e. warnings, errors, other key outputs) with a time stamp. This can make debugging easier, but can also help you pinpoint errors that do not throw obvious errors in the code (e.g. flawed image data that leads to unusual FE results).

Though there are many ways to accomplish this in Python, a utility function presents an excellent modular approach:

We can also initialize the time stamp after this point:

The registered FE method currently uses custom bone volumes from the XtremeCT II systems (insert bone image), so a material file makes it easy to identify the loading surfaces and assign material properties to the bone volume. For more details on creating material files, see Material definitions file.

Our material file will have the following format:

Now that we have a material file, we need to define a function that can read it for our model generator.

This function defines material names and identifiers based on the labels found in the bone volumes, then reads the material properties (i.e. modulus and Poisson’s ratio) from the material file. These properties will be assigned to the FE model via a material table in subsequent steps.

To read the file, the function searches through each line in the material file to identify and extract the elastic modulus and the Poisson’s ratio. After reading the file, the material properties of the bone and load surfaces are returned.

The most important function in this script is the model generator function, which follows this sequence:

We want to read a standard HR-pQCT image file (AIM), so we will use vtkboneAIMReader(). Note the use of the message() function we defined earlier.

Now we want to define the transformation matrix by which we will transform the our boundary conditions of our model. For longitudinal image data, this transformation will likely be in the form of a 4x4 affine matrix that is returned from 3D rigid-body image registration. This matrix contains the rotations and translations that are required to align two images (typically, follow-up images aligned to the baseline orientation). In the SCANCO implementation of image registration, the matrix is stored within a text file, as shown below.

To rotate the boundary conditions of any follow-up images (in a given image series) with respect to the baseline image, you will need to extract the 3x3 rotation matrix (the 9 elements in the upper-left of the 4x4 affine matrix) by adding the following lines to our CreateN88Model() function:

Now would be a good time to read the material file by calling the function we defined earlier. Remember to place this code segment within the CreateN88Model() function.

Since we have processed all our inputs, we can now create a finite element mesh from our image. We will apply a connectivity filter to our image, to make sure the model comprises a single object (i.e. with no floating pieces of bone or image noise in the model).

Let’s mesh our image:

To apply the material properties to the mesh we just created, we will need to make a material table. This process takes the elastic modulus and Poisson’s ratio that we read from the material file and assigns them to a linear-elastic material in the table.

We apply the material parameters to the mesh by compiling the mesh and material table into an initialized form of the finite element model, a “test base”.

Then we can tell Faim to identify the load surfaces, which are voxels/elements labeled with the number 127.

Applying transformed boundary conditions is relatively simple and requires only a few additional lines of code from a standard “non-registered” test.

First, we need to define a displacement vector for the model.

The displacement vector is typically oriented along the z-axis and has a magnitude of 0.102 mm, 1% of the image “thickness”, for a standard XtremeCT II bone model (to achieve 1% apparent strain).

Since we are transforming the boundary conditions of follow-up images, so they are aligned in the same direction as the baseline image, the displacement vector needs to be transformed with the rotation matrix we generated earlier.

Now we want to create a set of boundary conditions on the load surface, using the transformed displacement vector we just created. We will do this using the ApplyBoundaryCondition() function from vtkbone.

Before writing the n88model file, we want to define the model outputs that will be recorded when the model is solved. So, we will define the node sets from which we acquire our bone strength estimates (e.g. reaction force), specifically the top and bottom surfaces (loaded and fixed) of the model we just created.

Almost there! The last step in this function is to write the model to an n88model file. The filename will be defined by an input variable “+output_file+”.

The main function of this script processes the user’s inputs and sends them to the model generator function to create our n88model file. The input arguments will be defined with an argument parser, including mandatory and optional variables. We can also include (brief) help documentation for the user to clarify what input variables are required for this script to run properly.

The input arguments we will define are:

Written in code, this function will look like this:

To end the script, it is important to make sure all of these functions will actually be called by the program when it is run from the command line. This is done by a short code segment that checks to make sure the script has been run correctly from the command line.

Now we will put this method into practice. We will test this script, using a registered XtremeCT II image (distal radius) and performing a standard compression test to 1% apparent strain.

You will need the following files, provided in the tutorial files ( http://numerics88.com/downloads/ ):

The image file is an XtremeCT II image of a distal radius. The radius was imaged four times and each follow-up image (S2, S3, and S4) were misaligned from the baseline image (S1). The images were registered with 3D rigid-body image registration, identifying the the largest common volume (% overlap) between scans. For this radius, the volume was found to be 91.37% of the original bone volume. We will use this value to determine our compressive displacement for the test.

The transform file contains the affine transformation that is needed to transform the image (S4) to baseline (S1); therefore, each follow-up image will have the same volume (defined by the overlap) and different transformation files.

Let’s run our script in “help” mode (-h) to make sure we have all the inputs:

Here we can see the information that we added to the main function in our script with the argument parser. From this, we know that we are still missing the output_file and displacement arguments.

We will set the output filename to be the same as the input, just with the “n88model” extension instead of AIM, so it will look something like: REPO_004_RR_S4_REG_HOM_LS.n88model .

The volume of this image is 91.37% of the standard image volume (thickness 10.2mm), so we will scale the input displacement to be 91.37% of the standard displacement (0.102 mm): \(0.9137*0.102mm = 0.093mm\) . Calculating the displacement in this way provides a simplified “scaling factor” so that the apparent strain remains close to 1% for every scan in a study.

Now that we have all our inputs organized, let’s test our script:

If the script ran correctly, the output should look something like this:

A warning may appear while the script is run,

This is nothing to worry about and will not affect the results of the script.

You can now solve the model with faim:

Once your n88model is solved, the results will need to be processed using a slightly altered procedure than in a typical “non-registered” analysis. The non-registered analysis extracts the results as the vertical component only (i.e. the z-component, as in reaction force, displacement, etc.), since the test is conducted in exactly the vertical direction. With the registered approach, the results are transformed and so the test is no longer conducted in the vertical direction alone; therefore, the resulting estimates must include all components (i.e. x-, y-, z-components). To do so, we will take the magnitude of a result (u) as:

If the model used in our example above was solved correctly, the reaction force results (see REPO_004_RR_S4_REG_HOM_LS_analysis.txt) should look something like this:

Let’s say that we want to find the “registration-adjusted” reaction force estimate. Instead of using the total z-component of the reaction force (-8.214 kN), we want to calculate the total reaction force magnitude as:

This operation accounts for the transformed boundary conditions and allows us to report precise estimates of bone strength from longitudinal data. The easiest way to implement this adjustment is to calculate the magnitude of an estimate after the results have been parsed into a delimited file or spreadsheet for statistical analysis in Excel, SPSS, or R.

To calculate failure load with n88pistoia we will need to use a modified command than would be used in a standard analysis:

From our example, this command should then look like:

The full results set from this example can be parsed from the results files (*_analysis.txt or *_pistoia.txt) using the following commands:

This concludes the tutorial! You should now be able to transform displacement boundary conditions using vtkbone and inputs from 3D image registration. The example included in this tutorial is representative of the analysis that can be performed on a standard single-stack (10.2 mm) radius or tibia image from XtremeCT II.