Data Files for “Measure Logic” Problem

MeasureLogic.p3dat

; fname: MeasureLogic.p3dat
;
; Exercise the Mechanical Ball integration to the Measure Logic with PFC. 
; Porosity is computed using measurement spheres with random radii. 
; Results for different model resolutions are compared against closed 
; form solutions.
;
; Itasca Consulting Group, Inc.
;
; ============================================================================

fish def define_meas_parameters(nmeas, pcntg)
; ARGUMENTS: 
;  - nmeas : integer value that indicates the number of randomly-placed 
;            probing measurement circles that will be created
;            You may vary this parameter to see the effects on the final result
;  - pcntg : 'percentage' argument of the MEASURE function
; 
  global p = pcntg
  global n = nmeas
end


fish def compute_meas_center(r)
;Randomly computes the coordinates of the measurement sphere center
;ARGUMENTS:
;  - r : radius of the measurement sphere
;
  local xcmin = xmin+r
  local xcmax = xmax-r
  local ycmin = ymin+r
  local ycmax = ymax-r
  local zcmin = zmin+r
  local zcmax = zmax-r
  global xc = xcmin + (xcmax - xcmin)*math.random.uniform
  global yc = ycmin + (ycmax - ycmin)*math.random.uniform
  global zc = zcmin + (zcmax - zcmin)*math.random.uniform
end


fish define create_meas_regions
  measradius = array.create(n)
  local maxmeasradius = math.min((0.99*(zmax-zmin)),  ...
                        math.min((0.99*(xmax-xmin)),  ...
                        (0.99*(ymax-ymin))))
  loop ii(1,n)
    measradius(ii) = ...
          math.max(maxmeasradius/15,math.random.uniform*(0.5*(maxmeasradius)))
  end_loop
  loop ii(1,n)
    command
      @compute_meas_center([measradius(ii)])
      measure create position-x @xc position-y @yc position-z ...
                                @zc radius [measradius(ii)] tolerance @p
    end_command
  end_loop
end

fish define check_error(val,t,averageradius)
; Compares porosities computed in the measurement circles
; with analytical values
; ARGUMENTS: 
;   - val: exact value of the porosity
;   - t: name of the table with the final numerical results
;
  global result = table.create(t)
  local it = 0
  local averageerror = 0
  loop foreach local mp measure.list
    it = it + 1
    local res = measure.porosity(mp)
    local diff = 100.0 * (res - val) / val
    averageerror = averageerror + math.abs(diff)
    table(t,math.abs(diff))=measradius(it)/averageradius
  endloop
  global check_error = averageerror/it
end

; ============================================================================
; EoF: MeasureLogic.p3dat
Top

PolydisperseAssembly.p3dat

; fname: PolydisperseAssembly.p3dat
; Build a polydisperse granular assembly and compute porosity 
; using the measure logic in PFC
; ============================================================================
;

model new
call 'MeasureLogic.p3dat'

fish def poly_assembly(size)
; Create a polydisperse packing in a square domain
; ARGUMENTS:
;  - size : size of the domain
;
  global xmin       = 0.0
  global ymin       = 0.0
  global xmax       = size
  global ymax       = size
  global zmin       = 0.0
  global zmax       = size
  
  global width      = xmax-xmin  ; width of box
  global height     = ymax-ymin  ; height of box
  global depth      = zmax-zmin  ; depth of the box
  global tot_vol    = (width)*(height)*(depth)
  
  local num         = 10000   ; number of particles
  
  local rlo         = size/100
  local rhi         = 1.25*rlo
  local oversize    = rhi
  
  command
    model domain extent ([0.0-oversize],[width+oversize])  ...
                        ([0.0-oversize],[height+oversize]) ...
                        ([0.0-oversize],[depth+oversize])
    wall generate box (0.0, @width) (0.0, @height) (0.0, @depth) 
  end_command
 
  ;--- generate the balls and give them their properties
  command
    ball generate radius @rlo @rhi gauss box            ...
                         ([width*0.02], [width*0.98])   ...
                         ([height*0.02], [height*0.98]) ...
                         ([depth*0.02], [depth*0.98]) id=1,@num
  end_command
end

fish def averageradius
  local rad = 0
  loop foreach local bp ball.list
    rad = rad + ball.radius(bp)
  endloop
  global averageradius = rad/ball.num
end

fish def get_poros ;computes the exact value of porosity
  local sum = 0.0
  loop foreach local bp ball.list
    sum = sum + 1.33333* math.pi * ball.radius(bp)^3
  end_loop
  global get_poros = 1.0 - sum / tot_vol
end

@poly_assembly(5.0)

@get_poros
@define_meas_parameters(5000, 0.0)
@create_meas_regions

fish def display_error(porosity, s)
  global error = check_error(porosity,s,averageradius)
  local oo = io.out(string.build("Average Error = %1%)",error))
end
@display_error(@get_poros,'Polydisperse')

model save 'PolydisperseAssembly'

; ============================================================================
; EoF: PolydisperseAssembly.p3dat
Top

RegularAssembly.p3dat

; fname: RegularAssembly.p3dat
; Build a regular granular assembly and compute porosity 
; using the measure logic in PFC
; =============================================================================
;

model new
call 'MeasureLogic.p3dat'

fish define regular_assembly(size, nres,itype)
; Create a regular packing of mono-sized balls in a cubic domain
; ARGUMENTS:
;  - size : size of the domain
;  - nres : resolution of the packing
;  - itype: type of packing (0:simple cubic; #0: hexagonal)
; 
  global    xmin       = 0.0
  global    ymin       = 0.0
  global    xmax       = size
  global    ymax       = size
  global    zmin       = 0.0
  global    zmax       = size

  global    radius      = size / (2.0*nres)
  local     numballs    = nres^2

  command
    model domain extent @xmin @xmax @ymin @ymax @zmin @zmax
  end_command

  if itype = 0 then
    command
      ball generate rad @radius number @numballs cubic
    end_command
  else
    command
      ball generate rad @radius number @numballs hexagonal
    end_command
  endif
end

fish def display_error(porosity, s)
  global error = check_error(porosity,s,radius)
  local oo = io.out(string.build("Average Error = %1%)",error))
end

@define_meas_parameters(5000, 0.0) ;define measurement parameters

; Simple Cubic packing - low resolution
@regular_assembly(18.0, 9, 0)
@create_meas_regions
@display_error(0.4764,'cubic_lores')
model save 'cubic_lores'
;
; Simple Cubic packing - high resolution
ball delete
@regular_assembly(18.0, 30, 0)
@display_error(0.4764,'cubic_hires')
model save 'cubic_hires'
;
; Hexagonal packing - low resolution
ball delete
@regular_assembly(18.0, 9, 1)
@display_error(0.2595,'hexa_lores')
model save 'hexa_lores'
;
; Hexagonal packing - high resolution
ball delete
@regular_assembly(18.0, 30, 1)
@display_error(0.2595,'hexa_hires')
model save 'hexa_hires'

model save 'RegularAssembly'

; ============================================================================
; EoF: RegularAssembly.p3dat
Top

MeasureLogic.p2dat

; fname: MeasureLogic.p2dat
;
; Exercise the Mechanical Ball integration to the Measure Logic with PFC. 
; Porosity is computed using measurement circles with random radii. 
; Results for different model resolutions are compared against closed 
; form solutions.
;
; Itasca Consulting Group, Inc.
;
; ============================================================================

fish def define_meas_parameters(nmeas, pcntg)
; ARGUMENTS: 
;  - nmeas : integer value that indicates the number of randomly-placed 
;            probing measurement circles that will be created
;            You may vary this parameter to see the effects on the final result
;  - pcntg : 'percentage' argument of the MEASURE function
; 
  global p = pcntg
  global n = nmeas
end


fish def compute_meas_center(r)
; Randomly computes the coordinates of the measurement circle center
; ARGUMENTS:
;  - r : radius of the measurement sphere
;
  local xcmin = xmin+r
  local xcmax = xmax-r
  local ycmin = ymin+r
  local ycmax = ymax-r
  global xc = xcmin + (xcmax - xcmin)*math.random.uniform
  global yc = ycmin + (ycmax - ycmin)*math.random.uniform
end


fish define create_meas_regions
  global measradius = array.create(n)
  local maxmeasradius = math.min(0.99*(xmax-xmin),0.99*(ymax-ymin))
  loop ii(1,n)
    measradius[ii]= ...
        math.max(maxmeasradius/15,math.random.uniform*(0.5*(maxmeasradius)))
  end_loop
  loop ii(1,n)
    command
      @compute_meas_center([measradius[ii]])
      measure create id=@ii position-x @xc position-y @yc ...
                     radius [measradius[ii]] tolerance @p
    end_command
  end_loop  
end

fish define check_error(val,t,average_radius)
; Compares porosities computed in measurement circles with analytical values
; ARGUMENTS: 
;   - val: exact value of the porosity
;   - t: name of the table with the final numerical results
;
  global result = table.create(t)
  local it = 0
  local averageerror = 0
  loop foreach local mp measure.list
    it = it + 1
    local res = measure.porosity(mp)
    local diff = 100.0 * (res - val) / val
    averageerror = averageerror + math.abs(diff)
    table(t,math.abs(diff))=measradius[it]/average_radius
  endloop
  global check_error = averageerror/it
end

; ============================================================================
; EoF: MeasureLogic.p2dat
Top

PolydisperseAssembly.p2dat

; fname: PolydisperseAssembly.p2dat
; Build a polydisperse granular assembly and compute porosity 
; using the measure logic in PFC
; ============================================================================
;

model new
call 'MeasureLogic.p2dat'

fish def poly_assembly(size)
; Create a polydisperse packing in a square domain
; ARGUMENTS:
;  - size : size of the domain
;
  global xmin       = 0.0
  global ymin       = 0.0
  global xmax       = size
  global ymax       = size
  
  global width      = xmax-xmin  ; width of box
  global height     = ymax-ymin  ; height of box
  global tot_surface    = (width)*(height)
  
  local num        = 3500   ; number of particles
    
  local rlo        = size/200
  local rhi        = 1.25*rlo
  local oversize   = rhi

  command
    model domain extent ([0.0-oversize],[width+oversize]) ...
                        ([0.0-oversize],[height+oversize])
    wall generate box (0.0, @width) (0.0, @height)
  end_command
 
  ;--- generate the balls and give them their properties
  command
    ball generate radius @rlo @rhi gauss box ([width*0.01], [width*0.99]) ...
                         ([height*0.01], [height*0.99]) id=1,@num
  end_command
end

fish def averageradius
  local rad = 0
  loop foreach local bp ball.list
    rad = rad + ball.radius(bp)
  endloop
  global averageradius = rad/ball.num
end

fish def get_poros ;computes the exact value of porosity
  local sum = 0.0
  loop foreach local bp ball.list
    sum = sum + math.pi * ball.radius(bp)^2
  end_loop
  global get_poros = 1.0 - sum / tot_surface
end

@poly_assembly(20.0)

@get_poros
@define_meas_parameters(5000, 0.00)
@create_meas_regions

fish def display_error(porosity, s)
  global error = check_error(porosity,s,averageradius)
  local oo = io.out(string.build("Average Error = %1%)",error))
end
@display_error(@get_poros,'Polydisperse')

model save 'PolydisperseAssembly'

; ============================================================================
; EoF: PolydisperseAssembly.p2dat
Top

RegularAssembly.p2dat

; fname: RegularAssembly.p2dat
; Build a regular granular assembly and compute porosity 
; using the measure logic in PFC
; ============================================================================
;

model new
call 'MeasureLogic.p2dat'

fish define regular_assembly(size, nres,itype)
; Create a regular packing of mono-sized balls in a square domain
; ARGUMENTS:
;  - size : size of the domain
;  - nres : resolution of the packing
;  - itype: type of packing (0:simple cubic; #0: hexagonal)
; 
  global xmin       = 0.0
  global ymin       = 0.0
  global xmax       = size
  global ymax       = size

  global radius     = size / (2.0*nres) 
  local numballs    = nres^2

  command
    model domain extent @xmin @xmax @ymin @ymax
  end_command

  if itype = 0 then
    command
      ball generate rad @radius number @numballs cubic
    end_command
  else
    command
      ball generate rad @radius number @numballs hexagonal
    end_command
  endif
end

fish def display_error(porosity, tablename)
  global error = check_error(porosity,tablename,radius)
  local oo = io.out(string.build("Average Error = %1%)",error))
end

@define_meas_parameters(5000, 0.00) ;define measurement parameters

; Simple Cubic packing - low resolution
@regular_assembly(18.0, 9, 0)
@create_meas_regions
@display_error(0.2146,'cubic_lores')
model save 'cubic_lores.p2sav'
;
; Simple Cubic packing - high resolution
ball delete
@regular_assembly(18.0, 30, 0)
@display_error(0.2146,'cubic_hires')
model save 'cubic_hires.p2sav'
;
; Hexagonal packing - low resolution
ball delete
@regular_assembly(18.0, 9, 1)
@display_error(0.0931,'hexa_lores')
model save 'hexa_lores.p2sav'
;
; Hexagonal packing - high resolution
ball delete
@regular_assembly(18.0, 30, 1)
@display_error(0.0931,'hexa_hires')
model save 'hexa_hires.p2sav'

model save 'RegularAssembly.p2sav'

; ============================================================================
; EoF: RegularAssembly.p2dat
Top