Example 5j2: Augmented frictional compression of breast phantom (solving from intermediate soln)

This example demonstrates an augmented frictional contact compression on a breast phantom model. The model starts from gravity-loaded state. In order to cope with discontinuous nature of the compression, local C_0 continuity is enforced by adding versions at nodes where sudden change of shape is expected. The whole problem is divided into 4 load steps, and each load step is augmented 3 times to enforce the contact constraints. Solutions for each load step are saved as slave_soln_$i.ipinit ($i=1..4), and the Lagrangian multipliers of contact pressures are also saved as contact_pressure_$i.ipdata. However, this example illustrates how to run the simulation of the last step only, by reading in the solutions from the 3rd load step, as well as the contact pressure data file. The solution and contact pressure data files (slave_soln_test_4.ipinit, contact_pressure_test_4.ipdata) obtained by only solving for the last (4th) step can be compared with the solution/contact pressure data files (slave_soln_4.ipinit, contact_pressure_4.ipdata) obtained by solving from the beginning. As expected, there is no difference.


The comfile run by this example is as follows:

# 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 echo
$load_step=4;                      # set number of load step
$load_step_size=12/$load_step;     # set compression level at each load step
$BREAST=1;
$PLATE=2;
$START=4;     	      	      	  # specify which load step to solve.
      	      	      	      	  # set to 1 if solving from the beginning

fem def para;r;contact;example   # allocate memory
fem def coord;r;coord;example    # 3D space with C_0 continuity mapping option
fem def region;r;region;example  # define 2 regions
fem def bases;r;bases;example    # tricubic_hermite for geometry, trilinear for incompressibility
                                 # see .ipequa files

# define geometry for the breast phantom
# versions are added to allow local C_0 continuity

fem def node;r;symmetry_breast;example reg $BREAST
fem def elem;r;symmetry_breast;example reg $BREAST
fem def fibr;r;symmetry_breast;example reg $BREAST
fem def elem;r;symmetry_breast;example fib reg $BREAST

#fem export node;symmetry_breast as symmetry_breast reg $BREAST
#fem export elem;symmetry_breast as symmetry_breast reg $BREAST

# define geometry for the compression plate

fem def node;r;Plate;example reg $PLATE
fem def elem;r;Plate;example reg $PLATE
fem def fibr;r;Plate;example reg $PLATE
fem def elem;r;Plate;example fib reg $PLATE

#fem export node;plate as plate reg $PLATE
#fem export elem;plate as plate reg $PLATE

# Breast node, elem grouping #

fem group elem all external s3=1 as BACK_ELEM reg $BREAST
fem group elem all external s2=0 as SYMMETRIC_ELEM reg $BREAST
fem group elem all external s3=1 s2=0 as EDGE_ELEM reg $BREAST

fem group node xi3=1 external elem BACK_ELEM as BACK_NODE reg $BREAST
fem group node xi2=0 external elem SYMMETRIC_ELEM as SYMMETRIC_NODE reg $BREAST
fem group node xi2=0 xi3=1 external elem EDGE_ELEM as EDGE_NODE reg $BREAST

# Plate node, elem & face grouping #

fem group elem all external s3=0 as PLATE_BOTTOM_ELEM reg $PLATE
fem group elem all external s3=1 as PLATE_TOP_ELEM reg $PLATE

fem group node xi3=1 elem PLATE_TOP_ELEM external as PLATE_NODE reg $PLATE
fem group face all elem PLATE_BOTTOM_ELEM external xi3 low as PLATE_FACE reg $PLATE

# define 3D elasticity with contact, with neo-Hookean incompressible material

fem def equa;r;contact;example reg $BREAST
fem def mate;r;breast;example reg $BREAST

if($START==1){

  # READ IN GRAVITY SOLUTION TO START
  fem def init;r;GravitySolution;example reg $BREAST

} else {

  # READ IN previous solution and contact pressure values
  $temp=$START-1;
  fem def init;r;slave_soln_$temp;example
  fem def data;r;contact_pressure_$temp;example contact
} 

# COPY XP TO YP3
fem update solution geom YP_index 3 reg $BREAST

#####################################################################
# IMPORTANT : WHEN SOLVING FOR THE COTNACT PROBLEM
# INCLUDE THE EFFECT OF GRAVITY ONLY FOR THE 1ST LOAD STEP.
# OTHERWISE GRAVITY VECTOR WILL ACCUMULATE THOUGHOUT THE LOAD STEPS

# MAKE SURE IN THE PLATE IPINIT FILE, THE GRAVITY VECTOR SET TO 
# THE SAME AS THE GRAVITY VECTOR FOR THE BREAST AS GRAVITY VECTOR
# IS NOT REGION DEPENDENT
#####################################################################

# define 3D elasticity with contact, with neo-Hookean incompressible material

fem def equa;r;plate;example reg $PLATE
fem def mate;r;plate;example reg $PLATE
fem def init;r;plate;example reg $PLATE

# define contact parameter values
fem def contact;r;contact;example

# define mapping (only decouple d/dxi2 derivatives)
fem def mapping;r;mapping;example
    
# define LU solver

fem def solv;r;LU;example 

$ADD_GRAVITY=1;

# load steps
for ($i=1;$i<=$load_step;$i++){

  # Translate XP
  fem change nodes translate by 0,-$load_step_size,0 reg $PLATE

  if($i>=$START){

  
    $CONVERGED=0;
    while ($CONVERGED==0){ 
   
      print "                      \n"
      print " ==================== \n"
      print "    Load step = $i    \n"
      print " ==================== \n"
      print "                      \n"
    
      ######################################Projection
      # Update XP from YP
      fem up geom from sol to 1..3 reg $BREAST

      # Set contact Xi points on specified faces
      fem def xi;c contact_friction faces 5,10,33,37 points 4 boundary # 'boundary' option allows to put 
                                                                       # contact constraints along element
								       # boundaries
      
      # Define data at xi positions 
      fem def data;c from_xi

      # Update slave info
      fem update data field to slave

      # Project onto target face
      fem def xi;c closest_face faces PLATE_FACE reg $PLATE

      # Store projection gap in data fields
      fem update data field from gap reg $PLATE
      
      # Update master info
      fem update data field to master 

      # Place initial geometry YP(ny,3) into XP
      FEM up geom from sol YP_index 3 to 1..3 reg $BREAST

      #####################################Mechanics problem 
      
      # Solve finite elasticity/contact problem
      if ($ADD_GRAVITY==1){

        # for the 1st load step only, set inc=1 (default) to keep gravity force in effect 
        fem solve error 1.0D-07 inc 1 iterate 1000 
      } else {
      
        # for subsequent load steps, inc must be set to zero, otherwise the gravity vector will accumulate
        fem solve error 1.0D-07 inc 0 iterate 1000 # avoid accumulating of gravity vector by increment 0
      }

    }#CONVERGED

    # Reset gravity load  
    $ADD_GRAVITY=0;
    
#    fem export node;Slave_deformed_$i as Slave_deformed_$i field reg 1
#    fem export elem;Slave_deformed_$i as Slave_deformed_$i field reg 1
    
    # Save each load step solutions.
    if($START==1){    	# solution obtained by solving from the beginning 
    
      fem def init;w;slave_soln_$i
      fem def data;w;contact_pressure_$i contact
   
    } else {  	      	# solution obtained by solving from an intermediate load step
   
      fem def init;w;slave_soln_test_$i
      fem def data;w;contact_pressure_test_$i contact
   
    }
  }
}


Files used by this example are:

Name                            Modified     Size

example_5j2.com 14-Oct-2008 6.1k GravitySolution.exelem 14-May-2008 63k GravitySolution.exnode 14-May-2008 32k GravitySolution.ipinit 14-May-2008 191k LU.ipsolv 16-Aug-2010 2.2k LU.ipsolv.old 14-May-2008 2.0k Plate.ipelem 14-May-2008 536 Plate.ipelfb 14-May-2008 367 Plate.ipfibr 14-May-2008 2.1k Plate.ipnode 14-May-2008 14k bases.ipbase 14-May-2008 6.6k breast.ipmate 14-May-2008 2.2k contact.ipcont 25-Sep-2008 640 contact.ipequa 14-May-2008 2.1k contact.ippara 11-Oct-2008 5.8k contact_pressure_1.ipdata 11-Oct-2008 6.9k contact_pressure_2.ipdata 11-Oct-2008 6.9k contact_pressure_3.ipdata 11-Oct-2008 6.9k contact_pressure_4.ipdata 11-Oct-2008 6.9k contact_pressure_test_4.ipdata 11-Oct-2008 6.9k coord.ipcoor 11-Oct-2008 688 mapping.ipmap 11-Oct-2008 16k plate.exelem 14-May-2008 8.8k plate.exnode 14-May-2008 6.2k plate.ipequa 14-May-2008 2.1k plate.ipinit 14-May-2008 3.2k plate.ipmate 14-May-2008 2.2k region.ipregi 14-May-2008 93 slave_soln_1.ipinit 11-Oct-2008 191k slave_soln_2.ipinit 11-Oct-2008 191k slave_soln_3.ipinit 11-Oct-2008 191k slave_soln_4.ipinit 11-Oct-2008 191k slave_soln_test_4.ipinit 11-Oct-2008 191k symmetry_breast.exelem 14-May-2008 59k symmetry_breast.exnode 14-May-2008 31k symmetry_breast.ipelem 14-May-2008 6.9k symmetry_breast.ipelfb 14-May-2008 3.6k symmetry_breast.ipfibr 14-May-2008 7.1k symmetry_breast.ipnode 14-May-2008 74k

Download the entire example:

Name                      Modified     Size

examples_5_5j_5j2.tar.gz 17-Aug-2010 140k

Testing status by version:

StatusTestedReal time (s)
rs6000-aix
cm64SuccessWed Mar 4 01:33:02 2009105
cm64-debugSuccessTue Mar 3 03:32:53 20092281
x86_64-linux
cmSuccessSun Mar 6 00:03:10 201656
cm-debugSuccessSat Mar 5 00:14:19 2016383

Testing status by file:


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

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


CMISS Help / Examples / 5 / 5j / 5j2