Introductory Slides YouTube Video Tutorial
Contents
 Getting Started With Continuity
 Slide Show
Introduction to MultiScale Modeling in Continuity
Continuity is a problem solving environment for multiscale modeling in biomechanics, biotransport and electrophysiology. It is freely available for academic use under license from UCSD and downloadable at the Continuity download page. Several new updates are released per year. Although Continuity is generic and applicationneutral, we have developed it to support our research on cardiac physiology and pathophysiology, and therefore it is particularly well adapted to multiscale problems in cardiac electrophysiology and biomechanics including coupled, multiphysics, multiscale models of the heart. In this introductory guided tour, we focus on example applications in cardiac modeling.
Because multiscale biological models are most often valuable when used in close association with experimental investigation, Continuity incorporates extensive facilities for analyzing and fitting measured data, including biomedical images, and the ability to author and dynamically compile component models, such as constitutive equations for tissue properties, systems models for cellular dynamics e.g. ionic current models in electrophysiology, material coordinate definitions, and lumped parameter models of circulatory hemodynamics. A key objective and design priority of our continuing Continuity development efforts is that every aspect of a problem or simulation that you would consider to constitute the mathematical model or its biological parameters (as distinct from the numerical implementation) should be data that can be easily shared and reused. For that reason too, there is also a public Continuity model database that you can access directly from within the Continuity client. We are close to realizing that goal, though certain parts of the model, most notably the underlying partial differential equations (PDEs) themselves, are still currently expressed in Continuity source code.
At its core, Continuity is a PDE solver that implements finite element methods for nonlinear mechanics, reactiondiffusion and monodomain problems. A typical multiscale problem in Continuity spans three or more scales, typically comprising a system of ODEs at each point in the physical domain representing a network of biochemical or biophysical interactions, a constitutive law representing physical properties of the system that spatially couple these pointwise networks, a system of PDEs representing the conservation law governing the physics over the anatomy of the system, and a model of boundary conditions representing the interactions of the system with neighboring structures. For example, in a typical cardiac biomechanics problem, ODEs represent dynamic Markov models of myofilament activation and interactions in the cardiac myocyte, the constitutive equations describe the threedimensional anisotropic mechanical properties of the tissue with respect to local structural axes and how the myofilament interactions modify them dynamically, the governing PDEs for force and momentum balance are solved for an anatomic model of the heart subject to boundary conditions developed by a lumpedparameter model of the whole circulation. Similarly, a cardiac action potential propagation model comprises equations representing the dynamics of individual ion channels and whole cell ionic currents and potentials as systems of ODEs, that are spatially coupled via the monodomain PDEs by anisotropic conductivities and subject to boundary conditions representing the passive current spread through the torso.
Installing Continuity
Continuity can be downloaded for Windows, (Intel) Mac OS X and Linux platforms from the Continuity download page. Versions 7113 and later are preconfigured to allow the authoring and compilation of userdefined models such as material coordinate transformations and constitutive laws. Older versions require the Compiler Toolkit PlugIn for Windows or the SetUp Build Environment Script for Mac OS X found here. The entire history of releases since Fall 2005 is there in case you need to go back to an earlier working version in the event of a problem in a subsequent release. There is also a subversion (svn) source code repository for developers. Installation uses industrystandard mechanisms specific to each platform, such as dmg files for the MacOS and an msi point and click installer for Windows. See the help pages for detailed installation instructions. When you run Continuity, you will see the option to register yourself as a user from the splash screen. You don’t need to register to install or use Continuity, but access to the model database does require a registration key as it makes use of federally supported Internet resources at UCSD and the San Diego Supercomputer Center. We encourage you to register and provide us with information about your research interests so that we can improve the tools for your needs and provide useful statistical information to our federal sponsors.
The Continuity Architecture
Continuity consists of a client and a server. The client serves as the graphical user interface, and the server carries out the numerical calculations for simulation and problem solving. The installation process installs both on your computer. When you launch Continuity the graphical user interface (GUI) client starts and gives you the default option of running the client and server together on your local machine so that they behave as a single application. However, you can also use the Continuity client to connect to a remote server running on another host including a Linux cluster. There is no need for the client and server to be running on the same OS platform. The server has dynamically compilable and loadable numerical analysis functions for mesh handling, data fitting and image handling, biomechanics, electrophysiology, and transport processes. Correspondingly, the user can load clientside components corresponding to these separate modules that appear as new menus in the GUI client. Indeed, you will see the option to load certain modules directly in the splash screen. By default only the Mesh module is selected. You don’t have to decide at this stage which module you want to use. You can always load modules at any time. For a more complete summary of Continuity’s design and architecture, see the Continuity overview.
Using the Mesh Module to Make a Simple Ventricular Model
The Mesh module is the most commonly used module in Continuity, which is why it is active by default when you start Continuity. (Note that you can use the Preferences menu to change the default modules that are selected when Continuity is launched). One powerful and uncommon feature of Continuity is the ability to create and use meshes using curvilinear coordinate systems, such as cylindrical polar or spherical polar coordinates, as an alternative to rectangular Cartesian coordinates. You will soon be able to define your own curvilinear world coordinate system in Continuity, but for the meantime, there are four builtin world coordinate systems including Cartesian coordinates.
Prolate Spheroidal Approximations of Ventricular Geometry
Investigators working on ventricular mechanics in the 1960’s and 70’s, including Streeter, Hanna and Mirsky, approximated the left ventricle as nested ellipses of revolution. Here we show data from Streeter and Hanna (1973), who fitted truncated ellipses of revolution to measurements of left ventricular epicardial and endocardial major and minor axis dimensions in canine hearts, isolated and fixed under various loading conditions.
Just as using cylindrical coordinates simplifies the description of tubular objects, the natural coordinate system for ellipsoidal shapes is prolate spheroidal (elliptichyperbolicpolar) coordinates (λ,μ,θ). In prolate spheroidal coordinates, the longitudinal (μ) coordinate forms families of ellipses, each with a common focal point at x=d. The radial (λ) coordinates form hyperbolae. Streeter and Hanna showed that, the focal distance is remarkably constant both between the epicardial and endocardial surfaces and between systole and diastole in the dog left ventricle. Based on their data we will create a prolate spheroidal mesh in Continuity using a focal length of 3.7 cm. From within the Continuity Mesh menu, choose Edit→Coordinates…, select Prolate Spheroidal in the popup menu, enter 3.7 for the focus position and click OK to submit.

Continuity Command Summary
From within the Mesh menu, choose Edit→Coordinates…
Select Prolate Spheroidal in the popup menu
Enter 3.7 for the focus position
Click OK to submit.
From the major and minor radii reported by Streeter and Hanna, the focus, and the equations for prolate spheroidal coordinates we can easily compute the λ coordinates of the epicardial and endocardial surfaces, as shown in the table. Streeter and Hanna also defined a truncation factor as the fraction of the equatorapex distance between the base and the equator. Their estimate of nearly 0.5 translates to a μ coordinate at the base of 120°.
Finite Element Interpolation
We’ll use these values to specify nodes of our mesh, but first we need to choose a suitable finite element basis function to interpolate the geometric variables. Unlike finite difference methods, which discretize a continuum into a finite number of points and approximate the governing differential equations using differences, finite element methods approximate the continuum using a finite number of functions, usually piecewise polynomials spanning subdomains (elements) with parameters defined at the vertices (nodes). The governing equations are represented using a weak, integral form that gives rise to sums instead of differences. This facilitates the finite element assembly of the global system of equations by addition of the element matrices representing the discretization of the integral equations over each subdomain.
Continuity supports a variety of basis functions for isoparametric finite element interpolation. The simplest type of elements are piecewise linear Lagrange elements. Note that within each element, there is a local finite element (ξ) coordinate that always varies from 0 at the first node to 1 at the second. The parameters of the linear interpolation are the values (u_{1} and u_{2}) of the function at the nodes, and thus the basis functions (φ_{1} and φ_{2})are the corresponding linear functions that when multiplied by u_{1} and u_{2}, respectively, and added together result in the linear interpolation between u_{1} and u_{2}. Note also that the term isoparametric here refers to the fact that rather than interpolate a variable y as a function of geometric variables x, we interpolate both as functions of ξ, so that y(x) does not itself need to be a function. The major advantage of finite element basis functions such as linear Lagrange interpolation schemes is that, since the parameters of the interpolation are the value at global nodes shared by adjacent elements, (C_{0}) continuity of interpolated variables between neighboring elements is automatically ensured.
Simply by taking tensor products, we can use these simple 1D linear Lagrange basis functions to generate 2D quadralateral bilinear Lagrange elements as functions of ξ_{1} and ξ_{2} and 3D hexahedral trilinear Lagrange elements as functions of ξ_{1}, ξ_{2}, and ξ_{3}.
We are now ready to make our first mesh in Continuity, using the prolate spheroidal coordinates we calculated from Streeter and Hanna’s data. We will build a simple truncated, thickwalled prolate spheroidal mesh with a single linear 3D element. First, we need to select a trilinear Lagrange 3D basis function. Using the Edit→Basis… command from the Mesh menu, select Lagrange Basis Function→3D→LinearLinearLinear with 3 integration/collocation points for “Xi 1” (i.e. ξ_{1}), “Xi 2” (i.e. ξ_{2}), and “Xi 3” (i.e. ξ_{3}). Click Add and OK to submit.

Continuity Command Summary
Using the Edit→Basis… command from the Mesh menu, select Lagrange Basis Function→3D→LinearLinearLinear
Use 3 integration/collocation points for “Xi 1” (i.e. ξ_{1}), “Xi 2” (i.e. ξ_{2}), and “Xi 3” (i.e. ξ_{3}).
Click Add then OK to submit.
Nodes
Selecting Edit→Nodes… command from the Mesh menu brings up a large form where we can define the nodal values of geometric coordinates and other predefined variables such as fiber angles and fields. In our model, Coordinates 1, 2 and 3 are λ, μ and θ, respectively. Even though a trilinear element has eight element nodes we only need four nodes for this model, because it is axisymmetric, wrapping around on itself. We’ll see how to handle this when we define elements, shortly. Since we only selected one basis function, LinearLinearLinear Lagrange 3*3*3 probably already appears at the top of each coordinate column in the Nodes form. If not, select it from the menu for all three coordinates. Now, click the Insert Node button on the left three times to so that there are a total of 4 nodes. Then simply type in the nodal prolate spheroidal coordinates for each node that we computed for the canine left ventricle. Note that the default units for angular coordinates such as μ and θ is degrees, but that can be changed to radians using the radio button at the bottom of the form. We will also assign a label to each node, which will be useful later when we define boundary conditions. To do this enter a text string in the text entry field at the top of the form next to the Edit button to the right of the “Node Number:” display. Use the labels APEX_ENDO, BASE_ENDO, APEX_EPI, BASE_EPI for the nodes 14, respectively, as shown in the table. The underscore in the label serves the special purpose of enabling a compound label, which we will take advantage of later when we refine the mesh. The grayed out fields next to the unchecked boxes for each nodal parameter are not needed for linear elements. They represent derivative parameters that we will use later for higherorder Hermite interpolation schemes. Special note: When you are entering nodal labels, you must press return for the entry to be saved. The Import/Export/Graph button on the left panel provides a convenient way to import or export the nodes as tabdelimited tables that can be viewed or edited in the Continuity Table Manager or a spreadsheet program, such as Microsoft Excel. Click OK to submit the nodal coordinates and labels.

Continuity Command Summary
Using the Edit→Nodes… command from the Mesh menu
Make sure LinearLinearLinear Lagrange 3*3*3 appears in top of each coordinate
Insert 3 nodes, so there are a total of 4
Define each node with values and labels
Click OK to submit the nodal coordinates and labels
Elements
Next, we’ll use the Edit→Elements… command from the Mesh menu to create one element. By entering the eight node numbers for each trilinear element, we both define the elements and the directions of their local ξ coordinates. The graphic in the form helps remind us of the numbering convention in which nodes are entered starting at (ξ_{1},ξ_{2},ξ_{3})=(0,0,0) and ending with the node at (ξ_{1},ξ_{2},ξ_{3})=(1,1,1). By convention for meshes defined in curvilinear coordinates, we make ξ_{1} coincide (to the extent possible) with the polar (θ) axis, ξ_{2} with the “longitudinal” (z, φ or μ) axis, and ξ_{3} with the “radial” (r or λ) axis, while also trying to keep the ξ coordinates righthanded by convention. Note that this has the odd effect of requiring that ξ_{1} coincide with the negative θ axis in prolate spheroidal coordinates because θ is the third prolate spheroidal coordinate by convention, but it is the second coordinate in cylindrical or spherical polar coordinates. However, in this example, ξ_{1} wraps all the way around, so the order doesn’t matter. Hence following these conventions, the eight element nodes of our single trilinear element are defined by the four global nodes (14) entering the node sequence 1,1,2,2,3,3,4,4 in the element form for element 1, where we have incremented first circumferentially (ξ_{1}), then longitudinally (ξ_{2}) from apex to base, and finally radially (ξ_{3}) from endocardium to epicardium. Note that in this very simple case, our world curvilinear (θ_{j}) coordinates and our local finite element (ξ_{i}) coordinates coincide and both are righthanded systems, though their coordinate order is different. In general, though, we will not be dealing with such simple shapes so we can not expect the ξ coordinates to follow along world coordinate axes, though we do generally try to maintain rightedhandedness. You should also enter a label for Element 1 (e.g. “LV”) into the Element Label field at the top of the form (remember to press return). Click OK to submit the elements.

Continuity Command Summary
Using the Edit→Elements… command from the Mesh menu
Enter the node sequence 1,1,2,2,3,3,4,4 in the element form for element 1
Enter a label for Element 1 (e.g. “LV”) into the Element Label field at the top of the form
Click OK to submit the nodal coordinates and labels
File→Save→Model to save your model
Calculating Mesh
Now we are ready to use our new mesh. The button on the toolbar is a shortcut to the Mesh→Calculate Mesh… menu choice. Click OK in the dialog box. Behind the scenes, Continuity also carried out a File→Send to Server command, which updates the Continuity server with any client changes and can be performed using the toolbar button. Nothing visible happens yet, but now you can render the mesh as a wire frame or surfaces using the Mesh→Render→Elements… command or the toolbar button. There is also a Mesh→Render→Nodes… command. There is a wide variety of graphical rendering controls in the View menu. For example, the View→Show→OpenMesh command lists the rendered objects and their attributes. The entry field next to “select” in the Information tab provides a way to highlight specific elements or nodes. Select the nodes object from the list of graphical objects on the left to try this.

Continuity Command Summary
Calculate the mesh, using either or Mesh→Calculate Mesh…
Render elements using either or Mesh→Render→Elements…
Mesh Volumes
The Mesh→List→Elements… command will give you information on element volumes, faces, and lines. Check that the volume of the wall is being integrated accurately. If necessary, you can adjust the number of Gaussian quadrature integration points in the Edit→Basis… form to increase the accuracy of numerical integration, at the expense of more calculations per element. Although the left ventricular cavity itself is not part of the mesh, there are also ways to integrate the cavity volume bounded by defined surfaces in Continuity. Choose the Mesh→List→Volume… command. This command calculates volumes bounded by a surface surrounding a single origin. To get an accurate cavity volume, we need to place the origin at the plane of the endocardial base. We can use the information in the geometric data table to compute this position. Y and Z are zero but X is negative and equal to hxa, which is 0.45×4.2=1.89 for the endocardium at enddiastole. For (X,Y,Z) under “Define the Origin” enter (1.89,0.0,0.0). Use the pulldown menu next to the “Element List” field and you should see your previously defined element label LV. Select it. Choose LV Endocardium from the menu underneath “Wall”, then click the Insert Surfaces button. Then select Epicardium from the “Wall” menu and Insert Surfaces again. Click OK and the cavity and wall volumes will be displayed. Because the calculation is completely different, the wall volume will probably not be exactly the same as the mesh volume displayed previously, but it should be fairly close.

Continuity Command Summary
Mesh→List→Elements…
Mesh→List→Volume…
“Define the Origin” enter (1.89,0.0,0.0)
Select LV from “Element List” field
Select LV Endocardium from “Wall” field
Click the Insert Surfaces
Select Epicardium from the “Wall”
Click the Insert Surfaces again
Click OK
Ventricular Myofiber Architecture
The ventricular muscle fiber architecture has a characteristic pattern in most mammals that is surprisingly conserved. Roughly, the myofibers remain in planes almost parallel to the walls with helical pitch orientations varying from 60° to 90° counterclockwise from the circumferential axis (righthanded helices) in the subendocardium to 45° to 90° (left handed) in the epicardium. This familiar convention for fiber orientation is due to the quantitative measurements of Streeter and coworkers in the 1960’s, and for this fiber orientations reported using these angles are sometimes called Streeter angles. Most models assume that fibers lie parallel to the wall, though some investigators have measured small but significant outofplane transverse or imbrication angles, especially near the base and apex, which we will assume to be zero. It is important to remember that cardiac myocytes, unlike skeletal myocytes, are not really fibers at all; they are rodshaped cells with endtoend and sidetoside junctional connections and branches. Sections cut across the fibers reveal that the myofibers are organized by the perimysial collagen matrix into branching laminar myocardial sheets approximately four cells thick. Their orientations tend to have a greater local dispersion and commonly show bimodal distributions in some regions, but based on the predominant sheet orientations, it is now common in biomechanical and electrical models of the ventricles to use a local, orthonormal fibersheetnormal coordinate frame following the work of LeGrice and colleagues.
We will use published measurements of fiber and sheet angles to build the material coordinates in our model. These measurements show that muscle fiber and sheet orientations vary transmurally in the dog. A comprehensive analysis of fiber and sheet architecture in the canine ventricles was performed by LeGrice et al (1997). But even with our single element prolate spheroidal model, we can represent much of the variation in myofiber architecture. To do so we will need basis functions that are higher order than linear.
Quadratic Lagrange and Cubic Hermite Basis Functions
By placing a third node at the middle of each element (at ξ=0.5), we now have three parameters that we can use to interpolate second degree polynomials using the quadratic Lagrange basis functions. Note that each basis function is itself a quadratic having the property that it is one at the node corresponding to its own nodal parameter and zero at all of the other nodes. In the same way we can use four equally spaced nodes for cubic Lagrange interpolation, but a more useful strategy that requires only end (or vertex) nodes (at ξ=0 and 1) is to use the cubic Hermite basis functions. In 1D, these four basis functions correspond to two parameters at each of the two end nodes, one being the value of the interpolated variable and the second being its first derivative with respect to ξ.
A big advantage of cubic Hermite interpolation is not only that the elements can now support a higherorder variation of mesh variables, but that now, by sharing global nodal parameters, we can now achieve first derivative (C_{1}) continuity of mesh variables between elements. A practical limitation arises here though because the ξ coordinates themselves are not generally continuous. For example, if neighboring elements are different in length (which is a common and useful property of typical unstructured finite element meshes), then shared nodal derivatives with respect to ξ will not result in smooth slopes in Euclidean space. We correct this problem in Continuity by storing derivatives with respect to global arc lengths s at global nodes, and converting them into local ξ derivatives in each element using local element scaling factors, which are computed from the element arc lengths when you issue the Mesh→Calculate Mesh… command in Continuity.
In the same manner as we did for linear Lagrange basis functions, we can construct 2D and 3D quadratic Lagrange and cubic Hermite elements by taking tensor products. We can even combine different 1D basis functions in each direction, to create, e.g. a 6noded quadraticlinear Lagrange 2D element. Forming a 3D tricubic Hermite element results in an 8noded hexahedral element with 64 nodal parameters, comprising 8 derivative parameters per node. For the next step, we will need two new basis functions: a linearbicubic tensor product basis function with 4 degrees of freedom per node, that we will use to interpolate transverse angles, and a bilinearcubic basis function with two degrees of freedom per node for the sheet angles. Using the Edit→Basis… command from the Mesh menu, select Hermite Basis Function→3D→LinearLinearCubic with 3 integration/collocation points for each local Xi coordinate direction as before. Click Add and OK to submit.

Continuity Command Summary
Using the Edit→Basis… command from the Mesh menu, select Hermite Basis Function→3D→LinearLinearCubic
Use 3 integration/collocation points for each local Xi coordinate.
Click Add then OK to submit.
Fiber/Sheet Angles
For the purposes of demonstration here, we will approximate the measured transmural distributions of fiber and sheet angles in our 3D prolate spheroidal model by using the trilinear Lagrange basis functions that you created. Later we will learn about the data Fitting module of Continuity which enables us to fit nodal parameters of meshes directly to measurements of geometry or field variables. We will use a trilinear basis function to interpolate fiber angles and estimate nodal values from the plotted data. Since they are more variable than fibers we will only specify values for the sheet angles, which we will also approximate using trilinear Lagrange interpolation. Estimated parameters are tabulated here. To enter these values, let’s go back to the Edit→Nodes… form, which you select from the Mesh menu. Click on the Fiber Angles tab. Then select LinearLinearLinear 3*3*3 from the Select Basis Number dropdown menu below “Fiber Angle”, select LinearLinearLinear 3*3*3 under “Transverse Angle”, and LinearLinearLinear 3*3*3 under “Sheet Angle”. Now you will see that extra nodal parameter entry fields are checked. Use the table to enter the nodal value parameters (ignore s(3)).

Continuity Command Summary
Mesh→Edit→Nodes…
Click on the Fiber Angles tab
Select LinearLinearLinear 3*3*3 under Fiber Angle
Select LinearLinearLinear 3*3*3 under Transverse and Sheet Angle
Enter values
Click OK
Save model
Material Coordinates
Continuity uses four coordinate systems: Global Rectangular Cartesian coordinates (Y_{1},Y_{2},Y_{3}); curvilinear World coordinates (Θ_{1},Θ_{2},Θ_{31}); local finite element coordinates (ξ_{1},ξ_{2},ξ_{3}); and local material coordinates (x_{1},x_{2},x_{3}). The Mesh→Edit Material Coordinates command allows us to define an orthogonal transformation that rotates global Cartesian coordinates at any point in the mesh to the local orthonormal material coordinates. The standard predefined material coordinate transformation in Continuity uses the local ξ coordinates as intermediates to create a locally orthogonal body coordinate system with one axis aligned with ξ_{1}, one in the ξ_{1}ξ_{2} plane and the third perpendicular to the other two. These axes are then rotated through the fiber, transverse and sheet angles to generate material coordinates coinciding with the fiber, sheet and sheetnormal directions. The fiber, transverse and sheet angles are the local varying parameters of this model. They carry no intrinsic meaning to Continuity and are only labeled as such in the Nodes editor for historical reasons when this material coordinate transformation was builtin and hardcoded in Continuity. Within the context of Continuity, we call this definition of material axis transformation the Standard Material Coordinates model. The equations and parameters of this material model have already been defined and are stored in the Continuity database.
Go to Mesh>Edit>Material Coordinates…, if the model shown is not MatCoordStandard, hit ‘Cancel’ in that window. Then within the Continuity client, go to File→Library→Search… and enter MatCoordStandard in the top search window. Rightclick on the model that was found and select Load. It is IMPORTANT that you select the third option in the popup windows that starts with Retain current problem…. Then, if you go back to Mesh>Edit>Material Coordinates, you should now have the Standard Material Coordinate model equations and parameters. Compile the model under Compile…→Compile Code as: C Double Precision. When the compilation has completed, check the Parameters to make sure that the fiber_angle parameter is set to the Fiber Angle field. You can change the values of the parameters later without having to recompile the model. Click OK to exit the model editor.

Continuity Command Summary
Mesh→Edit→Material Coordinates…
Confirm that MatCoordStandard model is listed
If MatCoordStandard is not listed:
File→Library→Search…
Search for MatCoordStandard
Load model and select third option in Resolve Conflict form
Mesh>Edit>Material Coordinates…
Compile…→Compile Code as: C Double Precision
Check the Parameters to make sure that the fiber_angle parameter is set to the Fiber Angle field
Click OK
Save model
Solving a Simple Biomechanics Problem in Prolate Spheroidal Coordinates
Mesh Refinement
To solve a Biomechanics problem, we will first need to refine the mesh to allow for reasonably converged stress and strain solutions. Here we will refine the mesh into 5 elements longitudinally and 2 transmurally. Select Mesh→Refine… and under “New Element per old element in:” select 1 for “xi1”, 5 for “xi2” and 2 for “xi3”. Click OK. Then rerender the mesh after first remembering to use the Mesh→Calculate Mesh… command. Your mesh should now have ten elements. Note that the node numbers will have changed but the labels will be preserved. More importantly the new nodes will have inherited the labels in a logical way. The new node at the apex will have the label APEX, the new nodes at the base will have the label BASE, the new nodes along the inside surface will be labeled ENDO and the new nodes along the outside will be labeled EPI. Use the Edit→Nodes… command from the Mesh menu to open the Nodes form and confirm this.

Continuity Command Summary
Mesh→Refine…
Select 1 for “xi1”, 5 for “xi2” and 2 for “xi3” under “New Element per old element in:”
Click OK
Calculate Mesh
Render Elements
Save model
Biomechanics Problems
The Biomechanics module of Continuity uses a Newton iterative method to solve the continuum equations of motion for nonlinear, anisotropic, large deformation elasticity problems subject to traction and/or displacement boundary conditions. You can load the Biomechanics menu and module using the add module toolbar button, or from File>Load Modules. Various other special features of biomechanics problems, notably viscoelasticity, dynamic active force generation and boundary conditions generated by circulatory models can be included too. The necessary governing equations include the nonlinear kinematic straindisplacement relations for large deformations, conservation of linear and angular momentum equations from continuum mechanics, a nonlinear constitutive equation for hyperelastic materials, and in the case of incompressible models, the Lagrangian mass conservation law. The weighted residual formulation of the finite element equations uses the Galerkin finite element method, which is equivalent in this case to solving the virtual work integral equations of mechanics.
We will also change our basis functions for the undeformed and deformed geometric variables from trilinear to linearbicubic to allow a more continuous and higher order variation of stress and strain transmurally and longitudinally. The refine command will automatically generate the correct derivatives we need for the undeformed geometry. If you still have the Nodes form open, change the basis functions for “Coordinate 1”, “Coordinate 2” and “Coordinate 3” to LinearCubicCubic Hermite 3*3*3. Click OK. Observe that the nonzero parameters that the mesh refinement has calculated and inserted into the form coincide to the derivatives of the coordinates with respect to the directions in which they are changing, i.e. Coordinate 1 (Lambda) is changing in the xi3 direction, Coordinate 2 (mu) is changing wrt xi2.

Continuity Command Summary
Mesh→Edit→Basis… and select Hermite Basis Function→3D→LinearCubicCubic
Use 3 integration/collocation points for each local Xi coordinate.
Click Add then OK to submit.
Mesh→Edit→Nodes
Change Coordinates 1, 2, and 3 to be LinearCubicCubic Hermite 3*3*3
Click OK
Boundary Conditions
We will use the Biomechanics Edit→Boundary Conditions menu choice to assign boundary conditions. But first use the Biomechanics→Update→Initial conditions from undeformed nodes command to copy the undeformed nodal coordinates and basis functions to the initial conditions tab of the Boundary Conditions form. Next, we will use the Edit→Boundary Conditions form to assign nodal displacement boundary constraints. Click on the Deformed Coordinate 1 tab to specify constraints on the deformed Lambda coordinate, the Deformed Coordinate 2 tab to specify constraints on the deformed Mu coordinate, and the Deformed Coordinate 3 tab to specify constraints on the deformed Theta coordinate. For the node numbers, you can use the labels BASE.*, APEX.* and BASE_EPI to specify the list of nodes for each constraint as shown in the table. (note the “.*” is a wildcard option that will capture the nodes labeled BASE_EPI, APEX_ENDO etc).
Next we can assign pressure boundary conditions using the External Pressure tab of the Boundary Conditions Form. In the “Element Number:” field enter a commaseparated list of the element numbers on the inner surface of the model. These are elements 15. (You can use the tools in View→Show→OpenMesh to confirm this). Next to “Select Pressure Type”, select “Incremental” and assign the pressure increment to the inner faces. For passive filling of the left ventricle 0.1 kPa should work well.

Continuity Command Summary
Edit→Load Modules→Biomechanics
Biomechanics→Update and select Initial conditions from undeformed nodes
Biomechanics→Edit→Boundary Conditions
Click on Deformed Coordinate 1 to specify constraints on the deformed Lambda coordinate
Click on Deformed Coordinate 2 to specify constraints on the deformed Mu coordinate
Click on Deformed Coordinate 3 to specify constraints on the deformed Theta coordinate
Click on External Pressure and define inner surface elements
Select “Incremental” under “Select Pressure Type”
Specify pressure for inner faces
Click OK
Constitutive Model
To specify the material properties, we use the Edit→Constitutive Model command from the Biomechanics menu to open the Model Editor. Here you can edit an existing model or author and compile a new one provided the compilation environment is set up. One model that is in the Continuity database (BM_Ortho_of_Lagrangian_strains_comp_U8_sympy) is orthotropic exponential slightly incompressible strain energy function of Lagrangian strains proposed by Usyk and McCulloch for canine myocardium in which the strain energy is a function of strain components referred to local fibersheetnormal material axes. Of course, without recompiling we can change the parameters, of which there are seven, an overall scaling coefficient and one coefficient for each of the six strain components in fibersheetnormal coordinates (E_{ff}, E_{ss}, E_{nn}, E_{fs}, E_{fn}, E_{sn}) that appears in the exponent of this strainenergy function. However, we will start with a slightly simpler transversely isotropic exponential strain energy function in which the fiber direction of the tissue is stiffer than the transverse direction but there is no difference between the sheet and sheetnormal (crossfiber) directions. This can be found in the Continuity database by searching for: BM_TI_of_Lagrangian_strains_comp_sympy), and can be loaded by rightclicking and selecting load and the selecting the third option in the Resolve Conflict window. This model has only five parameters instead of eight for the orthotropic model. The last parameter is a “Bulk Modulus” of the “slightly incompressible” term in the model. It is relatively high and penalizes changes in volume as the tissue pressure rises. The higher this number, the less the volume of the tissue will change as the material deforms. A number that is too low can make the solution sensitive to the specific value. As the value gets higher the solution changes less but can be harder to obtain. We will start with the default transversely isotropic exponential model. Later we will switch to the orthotropic model. Leave the parameters in their default values. Click Submit.

Continuity Command Summary
File→Library→Search
Search for “BM_TI_of_Lagrangian_strains_compy_sympy”
Rightclick on found model and select “Load”
From Resolve Conflict window, select third option (“Retain current problem…”)
Biomechanics→Edit→Constitutive Model…
Verify model is loaded and click “OK” to exit
Check Model Variables Size Allocation
Before we solve the biomechanics problem, we should check that the dimensions of all the array variables in the model are correct to avoid simple computational errors, particularly indexing errors. Go to View→Edit Dimensions and ensure that All Dimensions are OK is displayed. If there is a message about variable sizes being too small or too large, click Apply Marked Recommendations to automatically change the allocated sizes of all incorrectly sized variables. If you click Show Details, you can view each variable and its current and recommended sizes.

Continuity Command Summary
View→Edit Dimensions
Verify all dimensions or OK, or if not, accept the recommended corrections by selecting “Apply Marked Recommendations”
Click “OK” to exit
Do a send
Calculate Mesh
Note:, if you do a send and get an error about surface_origin. Save your model, then restart Continuity and reload your model. If you do send after reloading your model, it will work fine.
Solving the Biomechanics Model and Displaying the Results
Now all that remains is to issue the Biomechanics→Solve command. If you specify ten time steps of 1.0 units each, then 10 previously specified pressure increments will be applied; the final pressure will be 1.0 kPa. Note that the solve form includes stopping criteria for the iterative solver based on the sum of the solution vector increments, the sum of the unconstrained residuals and the maximum number of iterations. For this problem, the default setting will work, but that is not always the case. If the solution is having difficulty converging, you can relax the stopping criteria by increasing the number of iterations and the tolerance for the sum of solution increments and sum of unconstrained residuals (usually up to 1.00e002). Decreasing the time step by orders of magnitude is also helpful.

Continuity Command Summary
Biomechanics→Solve Nonlinear
Change “Number of Steps” to 10
Click “Solve”
Note that the next command makes use of a Field Variable which must have a basis function assigned in the Nodes Form. For Field Variable 1 in the Field Vector 1 tab of the Edit→Nodes form, choose LinearLinearLinear Lagrange 3*3*3. Once the problem has solved we can use the Biomechanics→Render→Surface command: under “Variables” select T Cauchy Stress Tensor to visualize individual components of the stress or strain tensor solutions on specified element surfaces. If you want to change the color, go to View→Show Open Mesh …, Select Textured Field from the Rendering order list, then Click on Colors tab select RGB or BGR for Color Scheme:, then Close the window. If you specify ξ_{11}=0 or ξ_{11}=0.5 you will be able to see results rendered on longitudinal planes. Note that fiber strain (xx) is more uniform than the transverse (yy or zz) components.

Continuity Command Summary
Mesh→Edit→Nodes
On the Fiber Angles tab, make sure that LinearLinearLinear Lagrange 3*3*3 is listed under Fiber Angle
Click “OK”
Biomechanics→Render→Surface
select variable T Cauchy Stress Tensor
Click “OK”
Now save your model file, restart Continuity, reload the model and change the Constitutive model to the orthotropic exponential model. Change the Bulk Modulus in the parameter tab from 1000 to 400. For this model to converge, you will need to change some of the solve settings. Change the time step to 0.5 and take 20 steps instead of ten. Reduce Delta in the solution scheme tab from 1.e4 to 1.e7. Compare the solutions with the previous case.
Prebuilt models
Prebuilt models ready for the biomechanics solve can be found in the Continuity model database. To access them, go to File→Library→Search… and find GettingStartedBMtrans.cont6 for the transversely isotropic case and GettingStartedBMortho.cont6 for the orthotropic case. To load them, simply rightclick on the title of the model and select Load.
Fitting
Electrophysiology
Biomechanics
Slide Show
See a slide show with all the images from this article.