t   = 1 		! mm	thickness
kyy = 1/500		! 1/mm	curvature in the y direction
kxx = 1/2000		! 1/mm	curvature in the x direction
E   = 2.1e5		! N/mm2	Young's modulus
nu  = 0.35 		! - 	Poisson's ratio
P   = 250		! N	point load
d   = 5			! mm	diameter of the point load distribution area
lx  = 12*sqrt(t/kyy)	! mm	span in the x direction
*IF,kxx,EQ,0,THEN	! mm	span in the y direction
 ly = lx
*ELSE
 ly = 12*sqrt(t/ABS(kxx))
*ENDIF
nx  = 80		! -	number of elements in the x direction		
ny  = 80		! -	number of elements in the y direction
h   = d/20		! mm	element size in the middle
b   = h/4		! mm	load patch size (load is applied in small patches)

/PREP7
MPTEMP,,,,,,,, 		! material: isotropic
MPTEMP,1,0
MPDATA,EX,1,,E
MPDATA,PRXY,1,,nu
ET,1,SHELL281		! element type: 8 node quadrilateral
R,1,t,t,t,t, , , 	! element thickness

gx=lx/2-nx*h		! insert nodes
gy=ly/2-ny*h
mx=lx-(8*nx-6)*nx**2*h
my=ly-(8*ny-6)*ny**2*h
qx=4*nx*(1+3*nx-4*nx**2)
qy=4*ny*(1+3*ny-4*ny**2)
ty=1
*DO,j,0,2*ny
 ty=-ty
 tx=1
 *DO,i,0,2*nx
  tx=-tx
  *IF,tx+ty,LT,1,THEN
   x=i*(i*(3-2*i)*gx+mx)/qx
   y=j*(j*(3-2*j)*gy+my)/qy
   z=(x*x*kxx+y*y*kyy)/2
   N,,x,y,z,,,
  *ENDIF
 *ENDDO
*ENDDO

SHPP,OFF		! no warning aspect ratio
*DO,j,1,ny		! insert elements
 *DO,i,1,nx
  k1=1+(i-1)*2+(j-1)*(3*nx+2)
  k2=1+i+(2+(j-1)*3)*nx+(j-1)*2
  k3=1+(i-1)*2+j*(3*nx+2)
  E,k3,k3+2,k1+2,k1,k3+1,k2+1,k1+1,k2
 *ENDDO
*ENDDO

*DO,j,1,2*nx+1		! insert roller edges
 i=(3*nx+2)*ny+j
 NANG,i,1,0,NX(i)*kxx,,,,-NX(i)*kxx,-NY(i)*kyy,1
 D,i,,0,,,,UX,UY,,,,
*ENDDO
*DO,j,1,ny
 i=(j-1)*(3*nx+2)+2*nx+1
 NANG,i,,,,0,1,NY(i)*kyy,-NX(i)*kxx,-NY(i)*kyy,1
 D,i,,0,,,,UX,UY,,,,
 i=j*(3*nx+2)
 NANG,i,,,,0,1,NY(i)*kyy,-NX(i)*kxx,-NY(i)*kyy,1
 D,i,,0,,,,UX,UY,,,,
*ENDDO

*DO,i,1,2*nx+1		! insert symmetry boundary conditions
 D,i,,0,,,,UY,ROTX,ROTZ,,,
*ENDDO
*DO,j,1,ny+1
 D,(j-1)*(3*nx+2)+1,,0,,,,UX,ROTY,ROTZ,,,
*ENDDO
*DO,j,1,ny
 D,2*nx+2+(j-1)*(3*nx+2),,0,,,,UX,ROTY,ROTZ,,,
*ENDDO

sx=1			! insert point load
a=d/2-NX(2*sx+1)				! determine the range sx and sy of loaded elements
*DOWHILE,a
 sx=sx+1
 a=d/2-NX(2*sx+1)
*ENDDO
sy=1
a=d/2-NY((3*nx+2)*sy+1)
*DOWHILE,a
 sy=sy+1
 a=d/2-NY((3*nx+2)*sy+1)
*ENDDO

FCUM,ADD
n=0						! n = number of rings
a=d/2
*DOWHILE,a
 n=n+1
 a=a-b
*ENDDO
dr=d/2/n					! dr = ring width
*DO,k,1,n					! k = ring number
 m=0						! m = number of sectors
 a=3.1415/2*k*dr
 *DOWHILE,a
  m=m+1
  a=a-dr
 *ENDDO
 df=3.1415/2/m					! df = sector angle
 dP=P/(1/4*3.1415*d*d)*(2*k-1)/2*dr*dr*df	! dP = load on the sector
 ro=(12*k*k-12*k+4)/(6*k-3)*dr*SIN(df/2)/df	! ro = sector centre of gravity
 *DO,g,1,m					! g = sector number
  x=ro*cos((g-0.5)*df)
  y=ro*sin((g-0.5)*df)

  r=d*d						! find the closest node
  q1=1						! (look no further than sx or sy elements from the origin)
  *DO,j,1,sy
   *DO,i,1,sx
    q=(j-1)*(3*nx+2)+2*i-1
    a=(NX(q)-x)**2+(NY(q)-y)**2
    *IF,a,LT,r,THEN
     r=a
     q1=q
    *ENDIF
    q=(j-1)*(3*nx+2)+2*i
    a=(NX(q)-x)**2+(NY(q)-y)**2
    *IF,a,LT,r,THEN
     r=a
     q1=q
    *ENDIF
    q=(j-1)*(3*nx+2)+(2*nx+1)+i
    a=(NX(q)-x)**2+(NY(q)-y)**2
    *IF,a,LT,r,THEN
     r=a
     q1=q
    *ENDIF
   *ENDDO
  *ENDDO

  F,q1,FZ,dP					! move force to node and add
  F,q1,MX,dP*(y-NY(q1))
  F,q1,MY,dP*(NX(q1)-x)
*ENDDO
*ENDDO
FINISH

/SOLU				! compute
SOLVE 
FINISH

/POST1
*GET,w,NODE,1,U,Z		! deflection w under the point load
SHELL,TOP
*GET,sxxt,NODE,1,S,X		! stress sxx under point load in the top surface (z>0)
*GET,syyt,NODE,1,S,Y		! stress syy
*GET,sxyt,NODE,1,S,XY		! stress sxy
SHELL,BOT
*GET,sxxb,NODE,1,S,X		! stress sxx under point load in the bottom surface (z<0)
*GET,syyb,NODE,1,S,Y		! stress syy
*GET,sxyb,NODE,1,S,XY		! stress sxy
nxx=(sxxb+sxxt)*t/2		! membrane forces
nyy=(syyb+syyt)*t/2
nxy=(sxyb+sxyt)*t/2
mxx=(sxxb-sxxt)*t*t/12		! moments
myy=(syyb-syyt)*t*t/12
mxy=(sxyb-sxyt)*t*t/12

*CFOPEN,out,txt,,APPEND		! open file out.txt
*VWRITE,kxx,kyy,t,E,nu,P,d,lx,ly,nx,ny,h,w*E*t/P,mxx,myy,nxx,nyy
(F13.9,F13.9,F7.2,F10.0,F7.2,F10.0,F7.2,F10.0,F10.0,F5.0,F5.0,F6.2,F13.5,F13.5,F13.5,F13.5,F13.5)
*CFCLOS				! close out.txt
*UILIST,out.txt			! pop up out.txt

ETABLE,nxx,SMISC,1
ETABLE,nyy,SMISC,2
ETABLE,nxy,SMISC,3
ETABLE,mxx,SMISC,4
ETABLE,myy,SMISC,5
ETABLE,mxy,SMISC,6

FINISH