/* simple solid rocket motor */
options results

if ~show('ports', "AMI2D") then do
    address command 'run ami2d:bin/ami2d'
    do while ~show('ports', "AMI2D")
    end
end

address ami2d

call 'ami2d:rexx/axisym'

if ~show('l', "rexxmathlib.library") then do
    check = addlib('rexxmathlib.library',0,-30,0)
end

nn = 16
mn = 8
ne = nn/2 - 1
me = mn/2 - 1
b = 6.0
a = 1.0
h = 0.05
L = 10.0
E = 1500.0
nu = 0.4900
Ec = 10.0e6
nuc = 0.3000
P = -1000.0

'reset'
'iso(1,'E','nu')'
'iso(2,'Ec','nuc')'
g = 1.5
do i=0 to nn
    do j=0 to mn
        k = i*40 + j + 1
        r = a + (b-a)*pow(i/nn,g)
        z = -L*pow(j/mn,g);
        'node('k','r','z')'
    end
end
do i=1 to 2
    do j=0 to mn
        k = (nn + i)*40 + j + 1
        r = b + i*h/2
        z = -L*pow(j/mn,g);
        'node('k','r','z')'
    end
end
do i=0 to ne
    do j=0 to me
        k = i*20 + j + 1
        n = i*80 + j*2 + 1
        'quad8('k',1,'n','n + 2','n + 82','n + 80','n + 1','n + 42','n + 81','n + 40')'
    end
end
i = ne + 1
do j=0 to me
    k = i*20 + j + 1
    n = i*80 + j*2 + 1
    'quad8('k',2,'n','n + 2','n + 82','n + 80','n + 1','n + 42','n + 81','n + 40')'
end
do i=0 to (nn + 2)
    'dispz('i*40 + 1',0.0)'
end
do i=0 to me
    'press('i + 1',1,'P')'
end
j = me + 1
do i=0 to (ne + 1)
    'press('i*20 + j',2,'P')'
end
exit