Example 726: Material parameter estimation with an orthotropic Criscione material law through the CellMl environment for simple shear with forward solve

This example sets up a NF(xz)-shear for an eight element cube with two elements each direction to solve a material parameter estimation process with a forward solution procedure. It only validates the optmisation process. The material is assumed to be orthotropic and modelled by the Criscione material law through the CellMl environment, see also examples 5g* for further information about the implementation of an arbitrary constitutive law. The forward solution procedure is run to obtain synthetic force data of the top face of the cube which is later used for the objective function calculation, i.e. read in as the ipdata file. The optimiser then starts with different starting values and builds the gradient by calling the command file: "222_NF_Optimisation_Loop.com" which simulates the simple shear. As expected the optimiser converges to the values of theta51=0.8949 and theta53=0.2346 as stated in the ipcell file.

The comfile run by this example is as follows:

# The first part sets up the finite elasticity problem

fem define parameters;r;emech;example                   	# Reads in the parameter file
fem define coordinates 3,1					# Sets the coordinate system to be rectangular cartesian
fem define bases;r;TriLinLagr;example             		# Reads in a trilinear basis, only tri-and bilinear are 
								# needed since we use compressible material law
fem define;add bases;r;BiLinLagr;example 			# Reads in a bilinear basis, needed for the backend of CMISS
fem define node;r;222trilin_cube;example   			# Reads in the complete geometry including fibre orientation
fem define element;r;222trilin_cube;example        		
fem define fibre;r;222trilin_cube;example         
fem define element;r;222trilin_cube;example fibre    
fem group node 1,2,3,4,9,13,15,17,10 as BOTTOMFACE		# Groups the nodes of top and bottom face for the ipinit file
fem group node 5,6,7,8,11,14,16,18,12 as TOPFACE 
fem define equation;r;mechanics;example lock                    # Defines the finite elasticity equations
fem define material;r;CmissCellML;example			# Sets the connection to the CellMl environment
fem define initial;r;trilin_NearlyInc_222cube_SS_BC;example	# Defines the boundary conditions of a 0.5 shear in the NF-mode
								# Only the first three steps, i.e. 5%, 10% & 15% are used
								# for the optimisation procedure to speed up things.
fem define solve;r;trilin;example             			# Sets up a Newton Iteration with LU decompositio0n for 
								# the linearized problem
fem group elements 1..8 as ALL_ELEMENTS				# Groups All Elements, see below

# The Second part sets up the material law through the CellMl environment.
# For further information have look at examples 5g.

fem define grid;r;trilin_cube;example gauss class 2      	# Sets up grid points at Gauss points for the coupling
								# Make sure you use the same amount of Grid points 
								# as Gauss points in each direction
fem update grid geometry                      	
fem group grid element ALL_ELEMENTS as ALL_GRID_POINTS	   
fem define equation;r;cellml;example class 2                  	# This file sets up the proper conditions for the usage of
								# finite elasticity material laws
fem define material;d class 2

# These two files set up the nearly incompressible version of the orthotropic
# material law for laminar materials which is used here to simulate myocardium.
# The ipcell - file reads in the according .cml - file which is then compiled
# into fortran code by a CMISS library.
# This file reads in Green strains at Gauss points in the fibre coordinate system
# and ouputs 2nd Piola-Kirchhoff stresses.
# The material parameter values are NOT real parameters; they were obtained
# as a first rough estimate and shouldn't be used for real solutions.
# NOTE: theta1 should usually 100 to 1000 magnitudes higher than the other
# material parameters, if NEAR incompressibility is desired. 
# Otherwise use the compressible version of the Criscione law.

fem define cell;r;criscione_nearly_incompressible_legendre_two_parameters;example class 2      
fem define mate;r;criscione_nearly_incompressible_legendre_two_parameters;example cell class 2  

# The following commands are only set for the backend, they are not actually used
# during both the forward solve and the optimisation.

fem define initial;d class 2                        
fem define solve;r;cell-hmt-implicit-adams;example class 2 
fem solve class 2 to 0; #to num it is to 0 sec			# NOTE: It is important that the choice of sparsity versus 
								# non-sparisty in both ipsolv - files is consistent!

# The third part defines the actual optimisation procedure

#set output;SolveSolution3Steps on

# All data would be written into a file by the "set output;filename on" command.
# The following IF loop would run the forward problem and create the
# 222ForwardSolutionLoadPercent_$p - files for $solve=1. 
# Left out here for faster running. 

if ($solve) {
  for ($n=1;$n<=3;$n++)
        print "\n Increment: ", $n, "\n";
        fem solve iteration 10 increment 0.05 error 0.0000000000000001
        fem define initial;w;222ForwardSolutionLoadPercent_$p;example
  fem solve iteration 30 incr 0
  fem def init;w;222ForwardSolutionLoadPercent_3;example

# Reads in the "measured" data of the force applied to the TOPFACE
# with prescribed displacement during the shear experiment. In this case we only 
# read in data produced by the forward dolution

fem define data;r;NF_ExperimentalDataTest3;example

# The following loop writes the "measured" data into the proper
# array in CMISS. Make sure the number of residuals in the 
# ipopti - file match the number of residuals of the ipdata - file;
# 9 in our case.

for ($n=1;$n<=$number_load_steps;$n++)
	fem update residual data point $n index 1 residual 3*$n-2
	fem update residual data point $n index 2 residual 3*$n-1
	fem update residual data point $n index 3 residual 3*$n


# This >>symlink<< is set to current directory that full path doesn't need 
# to be specified in MaterialParameterOptimisation.ipopti
# For further information of symling have a look under
# >> perldoc -f symlink <<

symlink "$example/Optimisation_Loop.com","Optimisation_Loop.com"

# Defines the optimisation process. For our case we need to choose 
# Reaction differences for the objective function type. 
# The inquiry of the material parameter numbers is described at 
# the bottom of command file.

fem define optimise;r;MaterialParameterOptimisation;example

# Sets the output off.

# set output off

# The following commands were used to inquire the variable number
# of the material parameters within the CELL_RCQS_VALUE array 
# which are put into the MaterialParameterOptimisation.ipopti - file.
#fem inquire cell_variable theta51 return_variables THETA51ARRAY,THETA51IDX
#print $THETA51ARRAY 
#print $THETA51IDX
#fem inquire cell_variable theta53 return_variables THETA53ARRAY,THETA53IDX
#print $THETA53ARRAY 
#print $THETA53IDX
# These commands are used to watch the values of the above inquired
# variables at all Gauss points.
#fem list grid yqs index 26
#fem list grid yqs index 27

# For testing we compare the generated CMISS files
if ($TESTING) {
	fem def optimise;w;OptimisedMaterialParameters

Files used by this example are:

Name                                                            Modified     Size

example_726.com 22-May-2004 7.5k 222ForwardSolutionLoadPercent_10.ipinit 19-May-2004 27k 222ForwardSolutionLoadPercent_15.ipinit 19-May-2004 27k 222ForwardSolutionLoadPercent_5.ipinit 19-May-2004 27k 222trilin_cube.exelem 19-May-2004 12k 222trilin_cube.exnode 19-May-2004 5.4k 222trilin_cube.ipelem 19-May-2004 2.9k 222trilin_cube.ipelfb 19-May-2004 2.4k 222trilin_cube.ipfibr 19-May-2004 5.6k 222trilin_cube.ipnode 19-May-2004 6.4k BiLinLagr.ipbase 19-May-2004 1.1k CmissCellML.ipmate 19-May-2004 695 MaterialParameterOptimisation.ipopti 21-May-2004 1.4k NF_ExperimentalDataTest3.ipdata 19-May-2004 232 Optimisation_Loop.com 21-May-2004 1.8k OptimisedMaterialParameters.ipopti 21-May-2004 1.4k TriLinLagr.ipbase 19-May-2004 1.4k cell-hmt-implicit-adams.ipsolv 13-Apr-2007 2.6k cellml.ipequa 19-May-2004 1.5k criscione_nearly_incompressible_legendre_two_parameters.cml 19-May-2004 61k criscione_nearly_incompressible_legendre_two_parameters.ipcell 19-May-2004 8.5k criscione_nearly_incompressible_legendre_two_parameters.ipmatc 19-May-2004 2.1k emech.ippara 19-May-2004 5.9k mechanics.ipequa 19-May-2004 2.0k trilin.ipsolv 16-Aug-2010 2.2k trilin.ipsolv.old 13-Apr-2007 2.0k trilin_NearlyInc_222cube_SS_BC.ipinit 19-May-2004 1.4k trilin_cube.ipgrid 19-May-2004 610

Download the entire example:

Name                      Modified     Size

examples_7_72_726.tar.gz 27-Aug-2010 28k

Testing status by version:

StatusTestedReal time (s)
cmSuccessSun Mar 6 00:08:52 201612
cm-debugSuccessSat Mar 5 00:56:48 201644
cmFailureTue Aug 21 06:12:21 2007128
last breakTue Apr 17 02:12:00 2007127
last successSun Jun 6 01:55:00 200446
cm-debugFailureFri Aug 17 04:18:24 2007591
last breakTue Apr 17 02:27:00 2007578
last successFri Jun 11 01:38:00 2004204
cm-debug-clear-mallocFailureFri Aug 17 04:49:43 2007997
last breakSun Mar 19 03:20:00 20061172
last successSat Jun 5 01:57:00 2004237
cm-debug-clear-malloc7FailureFri Aug 10 05:09:07 20071087
last breakSun Mar 19 03:30:00 20061133
last successMon Jun 7 01:49:00 2004246
cm64FailureMon Aug 20 05:47:27 2007145
last breakTue Apr 17 02:58:00 2007130
cm64-debugFailureThu Aug 16 06:13:53 2007624
last breakTue Apr 17 03:16:00 2007613
cmFailureWed Mar 4 01:12:47 20091
last breakMon Aug 28 10:33:00 20061
cm-debugFailureWed Mar 4 01:13:47 20095
last breakThu Sep 13 02:09:00 20073
cm64FailureWed Mar 4 01:12:34 20091
last breakMon Aug 28 10:33:00 20061
cm64-debugFailureWed Mar 4 01:13:49 20093
last breakThu Sep 13 02:09:00 20073
cmSuccessSun Mar 6 00:01:22 20166
cm-debugSuccessSat Mar 5 00:04:47 201623

Testing status by file:

Html last generated: Sun Mar 6 05:50:31 2016

Input last modified: Thu Aug 26 14:19:53 2010

CMISS Help / Examples / 7 / 72 / 726