Example b361: Noble 98 - HMT 2D sheet with varying fibre orientation

This example uses the coupled Noble 98 - HMT cellular model within a 2D bicubic-Hermite sheet of tissue with the fibre angle varying from 10° at the lower edge of the sheet, to 90° in the middle of the sheet, and back to 10° at the upper edge (kind of visible in the images below) - but the mesh has been chopped down for testing purposes. The geometry shown below is available, just check the example command file for the file names. There is a 2:1 (fibre:cross fibre) conductivity ratio.

The model is electrically stimulated in the lower left corner. The sheet of tissue is free to deform as the propagating electrical wave triggers the active contraction at the cell level - with the boundary conditions that the extreme lower left node is held fixed in space and the extreme lower right node is free to slide along the horizontal axis.

The images below illustrate some results from this example, begining at time 0ms at the top left and running up to 550ms for the bottom right, in even steps of 50ms.

2D sheet simulation with varying fibre orientation

draw.com can be used to visualise the results of this simulation.

The comfile run by this example is as follows:

# Example b361:
# 2D sheet with N98-HMT used for active tension calculation
#

#
# This command file iterates through time in constant time steps - for
# each time step, the coupled cell model/distributed mechanics problem
# is iterated until a converged solution is obtained, such that the
# active tension generated at the cell level is in equilibrium with
# with the applied passive elastic loading at each point in time
#

# Look for the ##** sign to get the example running with the 8x8 element mesh.

# Set up some colour colour highlighting...
# To turn off all formating
$FBnone = "\033[0m";
# Foreground colurs
$Fblack   = "\033[30m";
$Fred     = "\033[31m";
$Fgreen   = "\033[32m";
$Fyellow  = "\033[33m";
$Fblue    = "\033[34m";
$Fmagneta = "\033[35m";
$Fcyan    = "\033[36m";
$Fwhite   = "\033[37m";
# Background colours
$Bblack   = "\033[40m";
$Bred     = "\033[41m";
$Bgreen   = "\033[42m";
$Byellow  = "\033[43m";
$Bblue    = "\033[44m";
$Bmagneta = "\033[45m";
$Bcyan    = "\033[46m";
$Bwhite   = "\033[47m";

# If the example path is not set, default to current directory
if (! defined $example) {
    $example = "./";
}
# Drop off the trailing / in the example path
$chopped = chop $example;
if ($chopped ne "/") {
    $example .= $chopped;
}

# Define the AndreSolve routine
if (! defined $andreSolveFile) {
    $andreSolveFile = $example."/andreSolve.com";
}
read command;$andreSolveFile;

#
# Main section ...
#

set echo;

$FORMAT = "ascii";

# for this example, we do not want to export data for every iteration
$output_iterations = 0;
# we do want to export cellular membrane potentials
$export_potentials = 1;
# we do not want to export any Gauss point information
$export_gauss = 0;
# and we do not want to list out any strains
$list_strains = 0;
# although, if we did we would want them at these Gauss points
@gauss_points = (1,38,39,42,43,22,23,26,27,64);

# Set the output directory
$output = ".";
if (! -d $output) {
    mkdir $output,0700;
}

# Create a log file for storing stuff
$LOGFILE = $output."/log.file";
CORE::open(LOG,">$LOGFILE") || die "Unable to open the log file ($LOGFILE): $!\n";
$date_string = time();
$date_string = localtime($date_string);
print LOG "Log file for Andre - created: $date_string\n\n";

# The two class definitions
$MECHANICS = 1;
$CELL = 2;

# Set to false to not solve the distributed mechanics
$SOLVE_DISTRIBUTED_MECHANICS = 1;

# Set to false to not solve the cell model
$SOLVE_CELL_MODEL = 1;

# Parameters for the distributed mechanics solution
$MECH_ERRTOL = 0.001;
$MAX_MECH_ITER = 20;

# The matrix file name to use for resetting the cell model
$MATRIXFILE = $output."/matrix";

# Set up the geometry 8x8 elements, bicubic Hermite basis
fem define parameter;r;$example/emech;
##**fem define node;r;$example/emech-8x8-cubic;
fem define node;r;$example/emech-4x4-cubic;
# Need a bi-cubic Hermite for geometry and fibres and a bi-linear
# for grid-based FEM
fem define base;r;$example/emech-cubic-linear;
##**fem define element;r;$example/emech-8x8-cubic;
fem define element;r;$example/emech-4x4-cubic;
# Fibres varying from 10 deg to 90 deg and back to 10 deg - with
# cubic Hermite interpolation
##**fem define fibre;r;$example/emech-8x8-109010;
fem define fibre;r;$example/emech-4x4-109010;
##**fem define element;r;$example/emech-8x8-cubic fibre;
fem define element;r;$example/emech-4x4-cubic fibre;

# Need to calculate the undeformed nodal derivatives
fem update node derivative 1 linear;
fem update node derivative 2 linear;
fem update node derivative 3 linear;

# Node/element groups to set BCs
# - the corners group will be fixed in space
# - the fixed_y group will be fixed in the y-direction
fem group node external as bdy_nodes
#fem group node 1,10,19,28,37,46,55,64,73 as corners
    fem group node 1 as corners
##**fem group node 9 as fixed_y
    fem group node 5 as fixed_y
##**fem group element 1..64 as ALL_ELEMENTS
    fem group element 1..16 as ALL_ELEMENTS

# Set up the finite elasticity
    fem define equation;r;$example/mech-cubic class $MECHANICS lock;
fem define material;r;$example/mech-cubic class $MECHANICS;
# use the cellular HMT active tension
fem define active;r;$example/mech-cubic class $MECHANICS;
fem define initial;r;$example/mech-cubic class $MECHANICS;
# Define the Newton iterative solver
fem define solve;r;$example/newton class $MECHANICS;

# Set up the N98-HMT cellular modelling - used to calculate the active 
# tension and wavefront propagation
fem define grid;r;$example/cell-18x18 class $CELL;
fem update grid geometry;
fem update grid metric;
#fem group grid xi1=0 oneoff as stimulus1_site;
fem group grid grid 1,2,3,10,11,12,19,20,21 as stimulus1_site;
##**fem group grid xi1 low element 5,23,21,29 as stimulus2_site;
fem group grid xi1 low element 3,11 as stimulus2_site;
fem group grid element ALL_ELEMENTS as ALL_GRID_POINTS;
# Use N98-HMT with grid-based FEM
fem define equation;r;$example/cell-n98_hmt class $CELL;
fem define material;r;$example/cell-n98_hmt class $CELL;
fem define cell;r;$example/cell-n98_hmt class $CELL;
fem define material;r;$example/cell-n98_hmt cell class $CELL;
fem update grid material class $CELL;
fem define initial;r;$example/cell-n98_hmt class $CELL;
# use the Adams-Moulton integrator
fem define solve;r;$example/cell-hmt-implicit-adams class $CELL;

# Export the undeformed geometry
fem export node;$output/sheet as sheet;
fem export element;$output/sheet as sheet;
# ... and define the deformed field ...
fem export node;$output/deformed_0 field as sheet;
fem export element;$output/deformed field as sheet;
# and the grid point numbers
fem export element;$output/numbers as sheet grid_numbers;

# Initialise the cell integration
fem solve class $CELL to 0;
# and write out the initial matrix file
fem write matrix;$MATRIXFILE matrices YQ,YQS $FORMAT;

# Tension index for grid (YQS) and Gauss (YG) variables
##fem inquire cell_variable Tension return_variables T_GRIDARRAY,T_GRIDIDX;
fem inquire cell_variable IsometricTension return_variables T_GRIDARRAY,T_GRIDIDX;
$T_GAUSSARRAY = "YG";
$T_GAUSSIDX = 1;
# !!!! These need to be greater than 6 so that when a update gauss
# !!!! stress/strain is done we don't loose anything
$T_FROM_CELLS_GAUSSIDX = 7;
$T_PREVIOUS_GAUSSIDX = 8;

# and the temporary storage of converged solutions in YP
$CONVERGED_SOLUTION_IDX = 11;

# Fibre extension ratio index for grid (YQS) and mechanics variables
fem inquire cell_variable ExtensionRatio return_variables L_GRIDARRAY,L_GRIDIDX;
$L_MECHIDX = 1;
# and previous fibre extension ratio index for grid (YQS) and mechanics variables
fem inquire cell_variable ExtensionRatioPrevious return_variables L_GRIDARRAY_PREV,L_GRIDIDX_PREV;
if (($L_GRIDARRAY ne $L_GRIDARRAY_PREV) || ($L_GRIDIDX == $L_GRIDIDX_PREV)) {
    die "Extension ratio current ($L_GRIDARRAY,$L_GRIDIDX) and previous ($L_GRIDARRAY_PREV,$L_GRIDIDX_PREV) arrays do not match!!";
}

# [Ca]i index for grid (YQS)
fem inquire cell_variable Cai return_variables CAI_GRIDARRAY,CAI_GRIDIDX;

# Vm index for grid (YQS)
fem inquire cell_variable Vm return_variables VM_GRIDARRAY,VM_GRIDIDX;

# Open the history files for tension, extension ratio, [Ca]i and Vm.
# One for the solutions for each iteration and one for the 
# coverged solutions
$ITERATIONS_UNIT = 21;
$CONVERGED_UNIT = 22;
fem open history;$output/cell-n98_hmt_iter write unit $ITERATIONS_UNIT variables yqs niqslist $T_GRIDIDX,$L_GRIDIDX,$CAI_GRIDIDX,$VM_GRIDIDX $FORMAT class $CELL;
fem open history;$output/cell-n98_hmt write unit $CONVERGED_UNIT variables yqs niqslist $T_GRIDIDX,$L_GRIDIDX,$CAI_GRIDIDX,$VM_GRIDIDX $FORMAT class $CELL;

# Initialise the distributed mechanics
# ?? Need to initialise YG ??
fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_FROM_CELLS_GAUSSIDX;
fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
fem update gauss gauss_field destination_field $T_PREVIOUS_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
if ($SOLVE_DISTRIBUTED_MECHANICS) {
    fem solve class $MECHANICS increment 0.0 iterate $MAX_MECH_ITER error $MECH_ERRTOL;
} else {
    $CONVERGED = 1;
}

#####################################################################
# NOTE:                                                             #
#   $CONVERGED is set in the CMISS routine NONLIN, either to 1 or 0 #
#     depending on whether convergence was achieved or not in the   #
#     non-linear solution.                                          #
#####################################################################


# save the intial converged solution
if ($CONVERGED) {
    fem update solution converged_reference save class $MECHANICS;
    if ($output_iterations && $export_gauss) {
        # Export the intial strains ??
        fem update gauss strain fibre components;
        fem export gauss;$output/gauss_strain_0 yg as sheet;
        # and stresses - first need to put the active tension back into
        # YG (the strains from the previous update have overwritten it)
        fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
        fem update gauss stress fibre components;
        fem export gauss;$output/gauss_stress_0 yg as sheet;
        # Need to replace the stress values that have been put into
        # YG(1..) with the initial active tension values
        fem update gauss gauss_field destination_field $T_GAUSSIDX source_field $T_FROM_CELLS_GAUSSIDX;
    }
} else {
    die "\n\nSomething is very, very wrong - see ya...\n\n";
}

# Set the initial fibre extension ratio for HMT grid problem
fem update grid extension_ratio component $L_MECHIDX $L_GRIDARRAY $L_GRIDIDX class $MECHANICS;

#
# The main time loop....
#
$tend = 10.0;
$tstep = 1.0;
$finish = 0;
$current_time = 0;
$previous_time = 0;
$history_time = 0;
set echo;
for ($time=0;($time<=$tend)&&(!$finish);$time+=$tstep)
{
    print LOG "*********************************************************************\n\n";
    print LOG "Time = $time\n\n";

    $previous_time = $current_time;
    $current_time = $time;

    # Perform the model solution for the current time step.
    # "sheet" is simply the graphical element group name to use for any
    # export to CMGUI files done during the solution process.
    ($converged,$num_iters) = &AndreSolve($current_time,$previous_time,$MATRIXFILE,"sheet");

    if ($converged) {
        print $Fcyan."Solve was successful, convergence was found in $num_iters iteration(s)$FBnone\n";
        print LOG "\n\nSolve was successful, convergence was found in $num_iters iteration(s)\n\n";
        # Write to the converged history file
        fem write history unit $CONVERGED_UNIT time $time variables yqs $FORMAT class $CELL;
        # Export the deformed nodes - need field on the export node to 
        # get the deformed geometry
        fem export node;$output/deformed_$time field as sheet;
        if ($export_gauss) {
            # Export the converged Gauss point strains
            fem update gauss strain fibre components;
            fem export gauss;$output/gauss_strain_$time yg as sheet;
            # Export the converged Gauss point stresses first need to 
            # put the active tension back into YG (the strains from the
            # previous update have overwritten it)
            fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_GAUSSIDX;
            fem update gauss stress fibre components;
            fem export gauss;$output/gauss_stress_$time yg as sheet;
            fem update gauss gridvars $T_GRIDARRAY $T_GRIDIDX $T_GAUSSARRAY $T_GAUSSIDX;
        }
        if ($list_strains) {
            # List out the strain at the specified Gauss points
            $gauss_string = "";
            foreach $gauss_point (@gauss_points) {
                fem list strain;$output/gauss_point_${gauss_point}_strain_$time gauss at $gauss_point fibre;
                $gauss_string .= "$gauss_point,";
            }
            chop $gauss_string;
            # and the grid point closest to them
            fem evaluate strain;$output/grid_point_strains_$time from cell_extension_ratios index $L_GRIDIDX at gauss $gauss_string nearest class $CELL;
        }
        if ($export_potentials) {
            fem export element;$output/potential_$time field as sheet class $CELL;
        }
    } else {
        print $Bred.$Fcyan."Solve was unsuccessful, convergence was not found with $num_iters iterations$FBnone\n";
        print LOG "\n\nSolve was unsuccessful, convergence was not found with $num_iters iterations\n\n";
        $finish = 1;
    }
    print LOG "=====================================================================\n\n";
} # end of time loop

# Close the cell history files
fem close history unit $ITERATIONS_UNIT $FORMAT class $CELL;
fem close history unit $CONVERGED_UNIT $FORMAT class $CELL;

# Compute the CMISS signal files from the history file
set echo;
# Vm
fem evaluate electrode;$output/Vm history $output/cell-n98_hmt from grid yqs iy $VM_GRIDIDX $FORMAT class $CELL;
fem evaluate electrode;$output/Vm_iter history $output/cell-n98_hmt_iter from grid yqs iy $VM_GRIDIDX $FORMAT class $CELL;
# [Ca]i
fem evaluate electrode;$output/cai history $output/cell-n98_hmt from grid yqs iy $CAI_GRIDIDX $FORMAT class $CELL;
fem evaluate electrode;$output/cai_iter history $output/cell-n98_hmt_iter from grid yqs iy $CAI_GRIDIDX $FORMAT class $CELL;
# tension
fem evaluate electrode;$output/tension history $output/cell-n98_hmt from grid yqs iy $T_GRIDIDX $FORMAT class $CELL;
fem evaluate electrode;$output/tension_iter history $output/cell-n98_hmt_iter from grid yqs iy $T_GRIDIDX $FORMAT class $CELL;
# and extension ratio
fem evaluate electrode;$output/extension_ratio history $output/cell-n98_hmt from grid yqs iy $L_GRIDIDX $FORMAT class $CELL;
fem evaluate electrode;$output/extension_ratio_iter history $output/cell-n98_hmt_iter from grid yqs iy $L_GRIDIDX $FORMAT class $CELL;

# Export to UnEMAP signal files

## Can only export from binary files, so don't do this unless you change the
## FORMAT variable to "binary".

# First the individual variables for all 729 grid points
##**fem define export;r;$example/cell;
#fem define export;r;$example/cell-4x4;
#fem export signal;$output/Vm electrode signal $output/Vm;
#fem export signal;$output/cai electrode signal $output/cai;
#fem export signal;$output/tension electrode signal $output/tension;
#fem export signal;$output/extension_ratio electrode signal $output/extension_ratio;
#fem export signal;$output/Vm_iter electrode signal $output/Vm_iter;
#fem export signal;$output/cai_iter electrode signal $output/cai_iter;
#fem export signal;$output/tension_iter electrode signal $output/tension_iter;
#fem export signal;$output/extension_ratio_iter electrode signal $output/extension_ratio_iter;

# Check for isovolumic contraction
fem list element total;
fem list element deformed total;

# close the log file
$date_string = time();
$date_string = localtime($date_string);
print LOG "\n\nLog file for Andre - closed: $date_string\n";
CORE::close(LOG);

quit

Files used by this example are:

Name                                Modified     Size


example_b361.com                    08-Oct-2002   14k
andreSolve.com                      08-Oct-2002   11k
cell-18x18.ipgrid                   15-Sep-2003  545 
cell-4x4.ipexpo                     12-Sep-2001  660 
cell-hmt-implicit-adams.ipsolv      13-Apr-2007  2.6k
cell-hmt-implicit-adams.ipsolv.old  12-Sep-2001  2.0k
cell-n98_hmt.ipcell                 12-Sep-2001   11k
cell-n98_hmt.ipequa                 26-May-2003  1.3k
cell-n98_hmt.ipinit                 12-Sep-2001  237 
cell-n98_hmt.ipmatc                 30-Jan-2003  673 
cell-n98_hmt.ipmate                 12-Sep-2001  1.7k
cell.ipexpo                         12-Sep-2001  660 
emech-4x4-109010.ipfibr             24-Sep-2001  7.1k
emech-4x4-cubic.ipelem              12-Sep-2001  5.1k
emech-4x4-cubic.ipelfb              12-Sep-2001  2.5k
emech-4x4-cubic.ipnode              12-Sep-2001   15k
emech-8x8-109010.ipfibr             24-Sep-2001  6.8k
emech-8x8-cubic.ipelem              12-Sep-2001   20k
emech-8x8-cubic.ipelfb              12-Sep-2001  9.7k
emech-8x8-cubic.ipnode              12-Sep-2001   49k
emech-cubic-linear.ipbase           12-Sep-2001  2.5k
emech.ippara                        21-Feb-2003  5.9k
mech-cubic.ipacti                   12-Sep-2001  253 
mech-cubic.ipequa                   02-May-2004  1.8k
mech-cubic.ipinit                   12-Dec-2002  1.6k
mech-cubic.ipmate                   12-Dec-2002  6.4k
newton.ipsolv                       16-Aug-2010  2.7k
newton.ipsolv.old                   13-Apr-2007  2.5k

Download the entire example:

Name                           Modified     Size


examples_b_b3_b36_b361.tar.gz  17-Aug-2010  330k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmSuccessSun Mar 6 00:12:51 201614
cm-debugSuccessSat Mar 5 00:39:53 201628
mips-irix
cmSuccessSun Aug 19 05:06:31 2007239
cm-debugSuccessWed Aug 15 04:09:08 2007410
cm-debug-clear-mallocFailureFri Aug 17 04:20:20 2007414
last breakTue Apr 17 06:08:00 2007415
last successSat Jun 25 01:40:00 2005207
cm-debug-clear-malloc7FailureTue Aug 21 05:38:02 2007415
last breakTue Apr 17 05:50:00 2007430
last successMon Jun 27 01:39:00 2005208
cm64SuccessSun Aug 19 05:10:50 2007226
cm64-debugSuccessTue Aug 21 03:26:42 2007421
cm64-debug-clear-mallocSuccessThu Apr 1 11:29:03 2004227
rs6000-aix
cmSuccessWed Mar 4 01:17:49 200916
cm-debugSuccessMon Mar 2 01:30:06 200983
cm64SuccessWed Mar 4 01:18:36 200919
cm64-debugSuccessTue Mar 3 01:33:58 200984
x86_64-linux
cmSuccessSun Mar 6 00:01:29 20168
cm-debugSuccessSat Mar 5 00:03:44 201615

Testing status by file:

Html last generated: Sun Mar 6 05:51:17 2016

Input last modified: Mon Aug 16 11:26:30 2010

CMISS Help / Examples / b / b3 / b36 / b361