Example 52c: Changing Dirichlet boundary conditions for quasi-static problems

This example demonstrates the use of the commands 'FEM change displacement_bc'. For quasi-static problems it is usually necessary to solve a sequence of problems just with a different set of boundary conditions. For example, in problems involving inverse dynamics, the movement of the bone is typically known from gait analysis (eg. during knee flexion or jaw movement). This movement also typically prescribes the movement of the muscle attachment area, and hence serves as boundary conditions for problms involving muscle mechanics. Instead of writing for each time step a seperate ipinit file, one can use this command to define the dirichlet boundary conditions within the comfile.

The comfile run by this example is as follows:


#
#	EXAMPLE 52c : CHANGING DIRICHLET BC FOR QUASI-STATIC PROBLEMS
#
#	This example demonstrates the use of the commands 'FEM change displacement_bc'.
#	For quasi-static problems it is usually necessary to solve a sequence of 
#	problems just with a different set of boundary conditions. For example, in 
# 	problems involving inverse dynamics, the movement of the bone is typically 
#	known from gait analysis (eg. during knee flexion or jaw movement). This movement 
#	also typically prescribes the movement of the muscle attachment area, and hence 
#	serves as boundary conditions for problms involving muscle mechanics. Instead
#	of writing for each time step a seperate ipinit file, one can use this command
#	to define the dirichlet boundary conditions within the comfile.
#


# 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;
}


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


FEM define para  ;r;  $example/MyRunParameter
FEM define coordinate 3,1
FEM define base  ;r;  $example/displacement

#
#	Read in data for a unit cube
#

FEM define node  ;r;  $example/UnitCube
FEM define elem  ;r;  $example/UnitCube
FEM define fiel  ;r;  $example/UnitCube
FEM define elem  ;r;  $example/UnitCube field
FEM define fibr  ;r;  $example/UnitCube
FEM define elem  ;r;  $example/UnitCube fibre

#
#	Simply defines a mechanics problem
#

FEM define equa  ;r;  $example/displacement 
FEM define mate  ;r;  $example/displacement

#
#	Writes/Exports reference configuration
#

FEM update solution from geometry
FEM update field from solution
FEM export elem;  $output/DefUnitCube0  field  element 1  as cube
FEM export nodes; $output/DefUnitCube0  field             as cube

#
#	Defines the initial condition. If using the commands 'FEM change disp', one does
#	specifies in the ipinit file only the nodes which are considered to be the
#	dirichlet boundary conditions (degrees of freedom which are allowed to change).
#	It's value is in this context not of importance, if one uses first the command
#	'FEM change disp init' and could therefor be any value. Typically, it is best 
#	to choose them as 0.
#

FEM define init  ;r;  $example/displacement


for ($i=1;$i<3;$i++) {

#
#	This command initialises all dirichlet boundary conditions to 0. This is important,
#	because the dirichlet bc always define the change from its current state to the next
#	one. This should be done at the beginning of each time step since otherwise the 
#	the changes to the dirichlet bc just get accumulated.
#
  FEM change disp init 2,4,6,8
 
  if ($i == 1) {

#
#	By 'FEM change disp rotate', one can specify te rotation angle, the axis and the nodes
#	which should be changed. Eg. the example below, sets the boundary conditions in way 
#	which is equivalent to rotating the right side of the reference cube by 22.5 degrees
#	around the x-axis. Note, only nodes which you specified as boundary nodes will be 
#	effected from the changes. 
#
#	By 'FEM change disp trans', one can specify a translation (add a vector) to the 
#	current/reference state. ' FEM change disp trans by  0.25,0,0' is the same as a 
#	strech of 0.25 units in the x-direction.
#	

      FEM change disp rotate     by  22.5           nodes 2,4,6,8  axis -1,0,0
      FEM change disp trans      by  0.25,0,0       nodes 2,4,6,8
  } else {

#
#	Instead of using explicit angles, one also can specify a 4x3 transformation matrix. The
#	matrix is defined in the following way:
#	          -                     -
#		  | r_11 r_12 r_13  t_x |
#	tmatrix = | r_21 r_22 r_23  t_y | 
#	          | r_31 r_32 r_33  t_z |
#	          -                     -
#	        = r_11,r_12,r_13,r_21,r_22,r_23,r_31,r_32,r_33,t_x,t_y,t_z
#
#	Hence, the command below is responsible for rotating the current state by -22.5 degrees
#	around the x-axis and moving it -0.25 units along the x-axis. The result is, the
#	rotation and translation from ($i==1) is reversed and we end up with a unit cube.
#

      FEM change disp tmatrix nodes 2,4,6,8  matrix  1,0,0,-0.25,0,0.92387953251129,-0.38268343236509,0,0,0.38268343236509,0.92387953251129,0	
  }    

#
#	Define the solver parameters
#

  FEM define solv  ;r;  $example/displacement
  FEM solve iter 20

#
#	Write/Export the current state for visualisation
#

  FEM update field from solution
  FEM export nodes;$output/DefUnitCube${i}  field as cube

}

#
#       The results can be visualised by the command View_Result.com
#       The time slider allows one to look at the different load steps.
#

Files used by this example are:

Name                     Modified     Size


example_52c.com          23-Aug-2006  4.6k
MyRunParameter.ippara    23-Aug-2006  6.0k
UnitCube.ipelem          23-Aug-2006  666 
UnitCube.ipelfb          23-Aug-2006  292 
UnitCube.ipelfd          23-Aug-2006  434 
UnitCube.ipfibr          23-Aug-2006  1.1k
UnitCube.ipfiel          23-Aug-2006   13k
UnitCube.ipnode          23-Aug-2006   12k
View_Results.com         23-Aug-2006  1.6k
displacement.ipbase      23-Aug-2006   12k
displacement.ipequa      23-Aug-2006  2.0k
displacement.ipinit      23-Aug-2006   49k
displacement.ipmate      23-Aug-2006  2.3k
displacement.ipsolv      16-Aug-2010  2.4k
displacement.ipsolv.old  13-Apr-2007  2.3k

Download the entire example:

Name                      Modified     Size


examples_5_52_52c.tar.gz  17-Aug-2010   16k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmSuccessSun Mar 6 00:03:37 20163
cm-debugSuccessSat Mar 5 00:27:20 201617
mips-irix
cmSuccessSun Aug 19 01:51:04 200720
cm-debugSuccessWed Aug 15 02:34:16 2007143
cm-debug-clear-mallocFailureMon Aug 20 06:14:16 2007138
last breakSat Aug 26 01:21:00 2006346
cm-debug-clear-malloc7FailureTue Aug 21 02:39:00 200712
last breakMon Aug 28 01:39:00 200642
cm64SuccessSun Aug 19 01:51:42 200720
cm64-debugSuccessTue Aug 21 02:38:46 2007161
rs6000-aix
cmSuccessWed Mar 4 01:09:36 20094
cm-debugSuccessMon Mar 2 01:22:19 200947
cm64SuccessWed Mar 4 01:09:36 20093
cm64-debugSuccessTue Mar 3 01:23:35 200945
x86_64-linux
cmSuccessSun Mar 6 00:01:04 20161
cm-debugFailureSun Mar 6 00:02:14 20167
last breakMon Dec 7 00:02:00 20157
last successSun Dec 6 00:03:00 20159

Testing status by file:

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

Input last modified: Mon Aug 16 11:19:14 2010

CMISS Help / Examples / 5 / 52 / 52c