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.
# 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&";