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 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.
$solve=0;
if ($solve) {
for ($n=1;$n<=3;$n++)
{
$p=5*$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.
$number_load_steps=3;
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
}
$ITERATION = 1;
# 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
optimise
# 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
}