TITLE 'Pore Interconnectivity'
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
COORDINATES
  CARTESIAN3('ro','Y','K')  
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
VARIABLES        
theta 
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
! SELECT     
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
DEFINITIONS    
eta=5

 L =  26 
 h = 0.20  

x=L*ro/2 
z=h*Y 

 P_in = 0.02 
 P_out = 0 
P=theta*(P_out-P_in)+P_in

k_z = 200 
k_x = K*k_z / ((2*h/L)^2)

blockage_start_bottom = -eta/L 
blockage_start_top = eta/L 

!Here's the integral that I'd like to solve
norm_flux = AREA_INTEGRAL(dY(theta), 'line')

V_x=-k_x*(2/L)*((P_out-P_in))*dro(theta)+P_in	
V_z=-k_z*(1/h)*((P_out-P_in))*dY(theta)+P_in	

! INITIAL VALUES
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

EQUATIONS        

theta: dro(dro(theta)) + (1/K)*dY(dY(theta))=0 

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
! CONSTRAINTS    { Integral constraints }
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
EXTRUSION K=0.00001,2

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
BOUNDARIES
  Region 'domain'
!mesh density cannot be increased very high for trial software
MESH_SPACING=.5

    Start 'outer' (-1,0)
      {walking around the sample starting at the bottom left corner which is unblocked }
    	value(theta)=1    line to (blockage_start_bottom,0)
		natural(theta)=0			line to (1,0)
		natural(theta)=0			line to (1,1) 
		value(theta)=0		line to (blockage_start_top, 1)
		natural(theta)=0			line to (-1,1)
		natural(theta)=0			line to close
Region 'line'
		Start 'outer1' (-1,0)
		line to (1,0)

!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
! TIME 0 TO 1    { if time dependent }
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
MONITORS         { show progress }
!---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
PLOTS            { save result displays }

!for eta by eta_o to eta_max
CONTOUR(theta) painted ON K=1

SURFACE(norm_flux) on Y=0

contour(norm_flux) ON Y=0

SUMMARY
Report ('fixed values in units of mm')
Report (eta) as 'Overlap'
Report (L) as 'Diameter'
Report (h) as 'Thickness'

END