! Clear data
/PREP7 
LSCLEAR,ALL 
CPDELE, 1   
CPDELE, 2   
CPDELE, 3   
CPDELE, 4   
CPDELE, 5   

/PREP7 

! Voltage couple top electrodes 
FLST,5,1408,1,ORDE,2
FITEM,5,1409
FITEM,5,-2816   
NSEL,S, , ,P51X 
CP, 1, VOLT, ALL
ALLSEL,ALL  
/REPLOT 

! Voltage couple bottom electrodes 
FLST,5,1408,1,ORDE,2
FITEM,5,1   
FITEM,5,-1408   
NSEL,S, , ,P51X 
CP, 2, VOLT, ALL
FINISH  
/SOL
!* 

! Static analysis 
ANTYPE,0

! Apply symmetry mechanical boundary conditions on XY = 0, YZ = 0 and ZX = 0 faces
FLST,2,3,5,ORDE,3   
FITEM,2,1   
FITEM,2,3   
FITEM,2,5   
DA,P51X,SYMM

FLST,2,1,5,ORDE,1   
FITEM,2,2   
/GO 
!*  

! Load in initial stress state
INISTATE, READ, file, ist


! Ground bottom electrodes   
ALLSEL,ALL  
/REPLOT 
FLST,2,1,1,ORDE,1   
FITEM,2,1   
/GO 
!*  
D,P51X,VOLT,0   

! Apply voltage on top electrodes
ALLSEL,ALL  
/REPLOT 
FLST,2,1,1,ORDE,1   
FITEM,2,1408
/GO 
!*  
D,P51X,VOLT,E_app 

! Convergence criteria
CNVTOL,F,,,,1E-6
CNVTOL,CHRG,,,,1E-6 
CNVTOL,VOLT,,,,1E-6 

		
            *DO,m,elem_1,elem_no ,1 
		        *IF,elem_mat,GT,mesh_size**3+2,THEN
                    *IF,(-2.969e-12*E_app**2+1.105e-5*E_app+9.003e1)*0.989,LT,RAND(1,100),THEN
                             MPCHG,elem_mat-max_elem-1,m
                    *ENDIF
				*ENDIF
			*ENDDO
			
			/STATUS,SOLU
			SOLVE 
				
			ESEL,S,MAT,,2,max_elem+1
            *GET,elem_unpoled_post,ELEM,0,COUNT
            ESEL,S,MAT,,n+1,meshdivs**3+n
            *GET,elem_poled_post,ELEM,0,COUNT
            ESEL,S,MAT,,1
            *GET,elem_air_post,ELEM,0,COUNT
            ALLSEL,ALL
            vf_air_post = elem_air_post
            vf_poled_post = elem_poled_post
            vf_unpoled_post = elem_unpoled_post
            con_res(icon,1) = vf_air_post
            con_res(icon,2) = vf_poled_post
            con_res(icon,3) = vf_unpoled_post
            *IF,icon,GE,2,THEN
            ic = icon - 1
            *VOPER,con_res(icon,4),con_res(ic,2),DIV,con_res(icon,2)
            *IF,con_res(icon,4),GT,0.95,AND,con_res(icon,4),LT,1.05,EXIT
            *IF,con_res(icon,4),EQ,1,EXIT
            *ENDIF





