Example 5j3: How to output contact parameters and visualise them in CMGUI

This example shows a simple frictionless contact problem between a soft cube under compression by a rigid (curved) plate. More importantly, it demonstrates how the contact parameters should be exported correctly (for debugging/visualisation purposes). Also, the command file "view.com" contains detailed explanation on how to display the exported files in CMGUI.


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

$CONT_PTS_1=4; 
$CONT_PTS_2=4; 
$load_step_size=2;
$BREAST=1;
$PLATE=2;

$output=0; # set this to 1 if you want to output/export files

fem def para;r;contact;example
fem def coord 3,1
fem def region;r;region;example
                                
fem def base;r;Base;example

fem def node;r;Cube;example
fem def elem;r;Cube;example
fem def fibr;r;Cube;example
fem def elem;r;Cube;example fibr

# Define field to store contact residual
fem def fiel;r;Cube;example
fem def elem;r;Cube;example field 

if($output){
  fem export node;Initial as Initial reg $BREAST
  fem export elem;Initial as Initial reg $BREAST
}

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

if($output){
  fem export node;Pad as Pad reg $PLATE
  fem export elem;Pad as Pad reg $PLATE
}

# node, elem & face grouping for B.C. #
fem group elem all external s1=1 as BACK_ELEM reg $BREAST
fem group node xi1=1 external elem BACK_ELEM as BACK_NODE reg $BREAST
fem group elem all external s1=0 as CONTACT_ELEM reg $BREAST
fem group face all xi1 low external elem CONTACT_ELEM as BREAST_FACE reg $BREAST
fem group face 18,76,128,165,197,206,215,220,227,234,241,244 as ONE
fem group face 7,52,65,110,121,148,156,192,43,101,143,185 as TWO

fem group elem 1001..1003 as MASTER_ELEM reg 2
fem group node 1001..1016 as MASTER_MOVE reg 2
fem group node xi3=1 elem MASTER_ELEM as MASTER_FIX reg $PLATE
fem group face all xi3 low external elem MASTER_ELEM as PLATE_FACE reg $PLATE

fem def equa;r;slave;example reg $BREAST
fem def mate;r;slave;example reg $BREAST
fem def init;r;slave;example reg $BREAST

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

fem def contact;r;contact;example

fem def solv;r;LU;example reg $BREAST

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

    #########################
    # Translate YP1 (current deformation)
    # Update XP from YP1
    fem up geom from sol to 1..3 reg $PLATE

    fem change nodes translate by $load_step_size,0,0 node MASTER_MOVE reg $PLATE
 
    # update YP1 from the translated XP
    fem update solution geometry substitute reg $PLATE
    #########################

    $j=1;
    $CONVERGED=0;
    while ($CONVERGED==0){ 
   
      ######################################Projection
      # Update XP from YP
      fem up geom from sol to 1..3 reg $BREAST,$PLATE

      # Set contact Xi points on specified faces
      fem def xi;c contact_point faces ONE,TWO points $CONT_PTS_1 by $CONT_PTS_2 boundary

      # Define data at xi positions 
      fem def data;c from_xi

      # Temporarily save the locations of the contact points on the slave face
      fem def data;w;junk

      # Update slave info
      fem update data field to slave

      # Project onto target face
      fem def xi;c closest_face search_start 1000 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 

      # export projection
      if($output){
      	$name=sprintf("projection.Load.%02d.Newton.%02d", $i, $j);
      	fem export data;$name as projection error offset 10000
      }
      
      # store actual locations of projected points on the master face
      fem def data;c from_xi

      # export normal and tangential vectors defined on the projected master face
      if($output){
      	$name=sprintf("projected.Load.%02d.Newton.%02d",$i, $j);
      	fem export data;$name as projected master offset 20000
      }
      
      # Place initial geometry YP(ny,3) into XP
      FEM up geom from sol YP_index 3 to 1..3 reg $BREAST,$PLATE

      #####################################Mechanics problem 

      # Solve finite elasticity/contact problem
      fem solve error 1.0D-08 iterate 25 reg $BREAST

      print " \n"
      print " ==================== \n"
      print "    Load step = $i \n"
      print " ==================== \n"
      print " \n"

      # Read-in temporary file to get correct geometric locations of contact points (on the slave)
      fem def data;r;junk;example

      # export contact forces as data file
      if($output){
        $name=sprintf("ContactPressure.Load.%02d.Newton.%02d", $i, $j);
        fem export data;$name as ContactPressure offset 30000 contact

        fem list node;ContactNodeForce.Load.$i.Newton.$j contact_forces
      }

      $j++; 

    }#CONVERGED

    if($output){

        ##### This is to visualise contact residual with current deformation #####
        ##### Must use XP instead of ZP #####
        # Update XP from YP
        FEM up geom from sol to 1..3 reg $BREAST

        # Write out contact residual to XP array 7 to 9
        FEM up geom from sol YP_index 6 to 7..9 reg $BREAST

        # Export final deformed meshes for current load step with contact residual
        fem export node;$name as Cube_deformed reg $BREAST
        fem export elem;$name as Cube_deformed reg $BREAST
    
        # Place initial geometry YP(ny,3) into XP
      	FEM up geom from sol YP_index 3 to 1..3 reg $BREAST
        ##### END #####

        $name=sprintf("plate_LoadStep_%02d",$i);
        fem export node;$name as Pad_LoadStep field reg $PLATE
      	fem export elem;$name as Pad_LoadStep field reg $PLATE
    }
}              

# write out contact residual, which is stored in YP(ny,6)
# WARNING: The contact residual written out here contains residual 
#          for reg 1 (under nc=1) and reg 2 (under nc=2) due to the 
#          different array setting for YP(ny,6) (no class is defined)
#          It is the same case for finite elasticity residual stored in
#          YP(ny,4), so reg 2 residual is stored under nc=2 in reg 1. 
fem list variable;contact_residual number 6

# save current solution
fem def init;w;conv_soln

## To see how to visualise these, look at view.com ##
#system "cmgui view.com&";


Files used by this example are:

Name                                      Modified     Size

example_5j3.com 17-Oct-2008 6.1k Base.ipbase 13-Oct-2008 6.6k ContactPressure.Load.04.Newton.09.exdata 13-Oct-2008 50k Cube.ipelem 13-Oct-2008 27k Cube.ipelfb 13-Oct-2008 18k Cube.ipelfd 13-Oct-2008 16k Cube.ipfibr 13-Oct-2008 24k Cube.ipfiel 13-Oct-2008 210k Cube.ipnode 13-Oct-2008 241k Initial.Frictionless.Load.04.exelem 13-Oct-2008 166k Initial.Frictionless.Load.04.exnode 13-Oct-2008 176k Initial.exelem 13-Oct-2008 166k Initial.exnode 13-Oct-2008 176k LU.ipsolv 16-Aug-2010 2.2k LU.ipsolv.old 13-Oct-2008 2.0k Pad.exelem 13-Oct-2008 13k Pad.exnode 13-Oct-2008 11k Pad.ipelem 13-Oct-2008 1.4k Pad.ipelfb 13-Oct-2008 944 Pad.ipfibr 13-Oct-2008 3.5k Pad.ipnode 13-Oct-2008 28k contact.ipcont 13-Oct-2008 626 contact.ippara 13-Oct-2008 5.8k conv_soln.ipinit 17-Oct-2008 658k junk.ipdata 13-Oct-2008 41k master.ipequa 13-Oct-2008 2.1k master.ipinit 13-Oct-2008 2.9k master.ipmate 13-Oct-2008 2.2k plate_LoadStep_05.exelem 13-Oct-2008 14k plate_LoadStep_05.exnode 13-Oct-2008 12k projected.Load.05.Newton.01.exdata 13-Oct-2008 68k projection.Load.05.Newton.01.exdata 13-Oct-2008 37k region.ipregi 13-Oct-2008 93 slave.ipequa 13-Oct-2008 2.1k slave.ipinit 13-Oct-2008 2.9k slave.ipmate 13-Oct-2008 2.2k view.com 13-Oct-2008 5.9k

Download the entire example:

Name                      Modified     Size

examples_5_5j_5j3.tar.gz 17-Aug-2010 233k

Testing status by version:

StatusTestedReal time (s)
rs6000-aix
cm64SuccessWed Mar 4 02:37:26 20091252
cm64-debugSuccessThu Nov 6 10:17:18 200819533
x86_64-linux
cmSuccessSun Mar 6 00:20:29 2016804
cm-debugSuccessFri Nov 7 04:14:51 20086724

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 / 5j3