*create,ceRVE

! This macro applies CE to a periodic hexaedral RVE
! The RVE must be centred at (0,0,0) and
! a node must exist at (0,0,0) 

a1=arg1
a2=arg2
a3=arg3

! APPLIED STRAINS
	eps11=arg4
	eps22=arg5
	eps33=arg6
!---Note: must enter epsij=1/2*gamma
	eps12=arg9/2
	eps13=arg8/2
	eps23=arg7/2  ! (x-y in model) 
!---Note: here recovers the other 1/2*gamma
	eps21 = eps12
	eps31 = eps13
	eps32 = eps23
! PERIODIC BOUNDARY CONDITION EQNS 2-DIR (X-DIR)
! DIR THEORY	ANSYS MODEL
!	1		Z
!	2		X
!	3		Y

!------------------------------------------------------
!   CREATE COMPONENT SET: periodic corners, edges and sides
!------------------------------------------------------

! select side +a1 and -a1 -----------------
nsel,s,loc,z,-a1-0.001,-a1+0.001
cm,a1n_nodes,node
nsel,s,loc,z,a1-0.001,a1+0.001
cm,a1p_nodes,node

! select side +a2 and -a2 -----------------
nsel,s,loc,x,-a2-0.001,-a2+0.001
cm,a2n_nodes,node
nsel,s,loc,x,a2-0.001,a2+0.001
cm,a2p_nodes,node

! select side +a3 and -a3 -----------------
nsel,s,loc,y,-a3-0.001,-a3+0.001
cm,a3n_nodes,node
nsel,s,loc,y,a3-0.001,a3+0.001
cm,a3p_nodes,node

! select edges a1 a2 ---------------------
cmsel,s,a1n_nodes
cmsel,r,a2n_nodes
cm,a1n_a2n_nodes,node
cmsel,s,a1n_nodes
cmsel,r,a2p_nodes
cm,a1n_a2p_nodes,node
cmsel,s,a1p_nodes
cmsel,r,a2n_nodes
cm,a1p_a2n_nodes,node
cmsel,s,a1p_nodes
cmsel,r,a2p_nodes
cm,a1p_a2p_nodes,node

! select edges a1 a3 ---------------------
cmsel,s,a1n_nodes
cmsel,r,a3n_nodes
cm,a1n_a3n_nodes,node
cmsel,s,a1n_nodes
cmsel,r,a3p_nodes
cm,a1n_a3p_nodes,node
cmsel,s,a1p_nodes
cmsel,r,a3n_nodes
cm,a1p_a3n_nodes,node
cmsel,s,a1p_nodes
cmsel,r,a3p_nodes
cm,a1p_a3p_nodes,node

! select edges a3 a2 ---------------------
cmsel,s,a3n_nodes
cmsel,r,a2n_nodes
cm,a3n_a2n_nodes,node
cmsel,s,a3n_nodes
cmsel,r,a2p_nodes
cm,a3n_a2p_nodes,node
cmsel,s,a3p_nodes
cmsel,r,a2n_nodes
cm,a3p_a2n_nodes,node
cmsel,s,a3p_nodes
cmsel,r,a2p_nodes
cm,a3p_a2p_nodes,node

cmsel,s,a1n_a2n_nodes
cmsel,a,a1n_a2p_nodes
cmsel,a,a1p_a2n_nodes
cmsel,a,a1p_a2p_nodes
cmsel,a,a1n_a3n_nodes
cmsel,a,a1n_a3p_nodes
cmsel,a,a1p_a3n_nodes
cmsel,a,a1p_a3p_nodes
cmsel,a,a3n_a2n_nodes
cmsel,a,a3n_a2p_nodes
cmsel,a,a3p_a2n_nodes
cmsel,a,a3p_a2p_nodes
cm,edges_nodes,node

! select corners a1 a2 a3 ---------------------
nsel,all
a1n_a2n_a3n_corner=node(-a2,-a3,-a1)
a1n_a2p_a3n_corner=node( a2,-a3,-a1)
a1n_a2n_a3p_corner=node(-a2, a3,-a1)
a1n_a2p_a3p_corner=node( a2, a3,-a1)
a1p_a2n_a3n_corner=node(-a2,-a3, a1)
a1p_a2p_a3n_corner=node( a2,-a3, a1)
a1p_a2n_a3p_corner=node(-a2, a3, a1)
a1p_a2p_a3p_corner=node( a2, a3, a1)
nsel,s,node,,a1n_a2n_a3n_corner
nsel,a,node,,a1n_a2n_a3p_corner
nsel,a,node,,a1n_a2p_a3n_corner
nsel,a,node,,a1n_a2p_a3p_corner
nsel,a,node,,a1p_a2n_a3n_corner
nsel,a,node,,a1p_a2n_a3p_corner
nsel,a,node,,a1p_a2p_a3n_corner
nsel,a,node,,a1p_a2p_a3p_corner
cm,corners_nodes,node
ALLSEL,ALL

NSEL,S,LOC,x,0
NSEL,R,LOC,Y,0
NSEL,R,LOC,Z,0
D,all,all
ALLSEL,ALL
CEDELE,ALL
ALLSEL,ALL
ceeq=0

! PERIODIC BOUNDARY CONDITION EQNS 2-DIR (X-DIR)
!
pos_node=
neg_node=
cmsel,s,a2n_nodes
cmsel,u,edges_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a2n_nodes
	cmsel,u,edges_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a2p_nodes
! 	get closest node from component neg_a2
	pos_node=node(-x_,y_,z_) 
	ceeq=ceeq+1
	ce,ceeq,(eps22*(a2*2)),neg_node,ux,-1,pos_node,ux,1	! x->2
	ceeq=ceeq+1
	ce,ceeq,(eps32*(a2*2)),neg_node,uy,-1,pos_node,uy,1	! y->3	
	ceeq=ceeq+1
	ce,ceeq,(eps12*(a2*2)),neg_node,uz,-1,pos_node,uz,1	! z->1
*enddo

! PERIODIC BOUNDARY CONDITION EQNS 3-DIR (Y-DIR)
!

pos_node=
neg_node=
cmsel,s,a3n_nodes
cmsel,u,edges_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a3n_nodes
	cmsel,u,edges_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a3p_nodes
! 	get closest node from component neg_a2
	pos_node=node(x_,-y_,z_) 
	ceeq=ceeq+1
	ce,ceeq,(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1	! x->2
	ceeq=ceeq+1
	ce,ceeq,(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3  
	ceeq=ceeq+1
	ce,ceeq,(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1
*enddo

! PERIODIC BOUNDARY CONDITION EQNS 1-DIR (Z-DIR)
!
pos_node=
neg_node=
cmsel,s,a1n_nodes
cmsel,u,edges_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a1n_nodes
	cmsel,u,edges_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a1p_nodes
! 	get closest node from component neg_a2
	pos_node=node(x_,y_,-z_) 
	ceeq=ceeq+1
	ce,ceeq,(eps21*(a1*2)),neg_node,ux,-1,pos_node,ux,1	! x->2
	ceeq=ceeq+1
	ce,ceeq,(eps31*(a1*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
	ceeq=ceeq+1
	ce,ceeq,(eps11*(a1*2)),neg_node,uz,-1,pos_node,uz,1	! z->1 
*enddo

! PERIODIC BOUNDARY CONDITION EQNS 23-DIR (XY-DIR)
!
pos_node=
neg_node=
cmsel,s,a3n_a2n_nodes
cmsel,u,corners_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a3n_a2n_nodes
	cmsel,u,corners_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a3p_a2p_nodes
! 	get closest node from component neg_a2
	pos_node=node(-x_,-y_,z_) 
	ceeq=ceeq+1
	ce,ceeq,(eps22*(a2*2))+(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1	! x->2   
	ceeq=ceeq+1
	ce,ceeq,(eps32*(a2*2))+(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
	ceeq=ceeq+1
	ce,ceeq,(eps12*(a2*2))+(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1
*enddo
!
pos_node=
neg_node=
cmsel,s,a3n_a2p_nodes
cmsel,u,corners_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a3n_a2p_nodes
	cmsel,u,corners_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a3p_a2n_nodes
! 	get closest node from component neg_a2
	pos_node=node(-x_,-y_,z_) 
	ceeq=ceeq+1
	ce,ceeq,(-eps22*(a2*2))+(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1	! x->2   
	ceeq=ceeq+1
	ce,ceeq,(-eps32*(a2*2))+(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
	ceeq=ceeq+1
	ce,ceeq,(-eps12*(a2*2))+(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1
*enddo

! PERIODIC BOUNDARY CONDITION EQNS 31-DIR (YZ-DIR)
!
pos_node=
neg_node=
cmsel,s,a1n_a3n_nodes
cmsel,u,corners_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a1n_a3n_nodes
	cmsel,u,corners_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a1p_a3p_nodes
! 	get closest node from component neg_a2
	pos_node=node(x_,-y_,-z_) 
	ceeq=ceeq+1
	ce,ceeq,(eps21*(a1*2))+(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1	! x->2   
	ceeq=ceeq+1
	ce,ceeq,(eps31*(a1*2))+(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
	ceeq=ceeq+1
	ce,ceeq,(eps11*(a1*2))+(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1
*enddo
!
pos_node=
neg_node=
cmsel,s,a1n_a3p_nodes
cmsel,u,corners_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a1n_a3p_nodes
	cmsel,u,corners_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a1p_a3n_nodes
! 	get closest node from component neg_a2
	pos_node=node(x_,-y_,-z_) 
	ceeq=ceeq+1
	ce,ceeq,(eps21*(a1*2))-(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1 	! x->2  
	ceeq=ceeq+1
	ce,ceeq,(eps31*(a1*2))-(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
	ceeq=ceeq+1
	ce,ceeq,(eps11*(a1*2))-(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1
*enddo

! PERIODIC BOUNDARY CONDITION EQNS 12-DIR (YZ-DIR)
!
pos_node=
neg_node=
cmsel,s,a1n_a2n_nodes
cmsel,u,corners_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a1n_a2n_nodes
	cmsel,u,corners_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a1p_a2p_nodes
! 	get closest node from component neg_a2
	pos_node=node(-x_,y_,-z_) 
	ceeq=ceeq+1
	ce,ceeq,(eps21*(a1*2))+(eps22*(a2*2)),neg_node,ux,-1,pos_node,ux,1	! x->2   
	ceeq=ceeq+1
	ce,ceeq,(eps31*(a1*2))+(eps32*(a2*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
	ceeq=ceeq+1
	ce,ceeq,(eps11*(a1*2))+(eps12*(a2*2)),neg_node,uz,-1,pos_node,uz,1	! z->1
*enddo
!
pos_node=
neg_node=
cmsel,s,a1n_a2p_nodes
cmsel,u,corners_nodes
*get,num_nodes,node,0,count,max
*do,i,1,num_nodes,1
	cmsel,s,a1n_a2p_nodes
	cmsel,u,corners_nodes
! 	neg_node is undefined; use lowest active node number when i=1
	*if,i,ne,1,then
		neg_node=ndnext(neg_node)
	*else
		*get,neg_node,node,0,num,min
	*endif	
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
	cmsel,s,a1p_a2n_nodes
! 	get closest node from component neg_a2
	pos_node=node(-x_,y_,-z_) 
	ceeq=ceeq+1
	ce,ceeq,(eps21*(a1*2))-(eps22*(a2*2)),neg_node,ux,-1,pos_node,ux,1 	! x->2  
	ceeq=ceeq+1
	ce,ceeq,(eps31*(a1*2))-(eps32*(a2*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
	ceeq=ceeq+1
	ce,ceeq,(eps11*(a1*2))-(eps12*(a2*2)),neg_node,uz,-1,pos_node,uz,1	! z->1
*enddo

! PERIODIC BOUNDARY CONDITION Corners
!
nsel,s,node,,a1n_a2n_a3n_corner
*get,neg_node,node,0,num,min
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
nsel,a,node,,a1p_a2p_a3p_corner
pos_node=node(-x_,-y_,-z_)
ceeq=ceeq+1
ce,ceeq,(eps21*(a1*2))+(eps22*(a2*2))+(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1	! x->2   
ceeq=ceeq+1
ce,ceeq,(eps31*(a1*2))+(eps32*(a2*2))+(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
ceeq=ceeq+1
ce,ceeq,(eps11*(a1*2))+(eps12*(a2*2))+(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1

nsel,s,node,,a1n_a2n_a3p_corner
*get,neg_node,node,0,num,min
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
nsel,a,node,,a1p_a2p_a3n_corner
pos_node=node(-x_,-y_,-z_)
ceeq=ceeq+1
ce,ceeq,(eps21*(a1*2))+(eps22*(a2*2))-(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1 	! x->2  
ceeq=ceeq+1
ce,ceeq,(eps31*(a1*2))+(eps32*(a2*2))-(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
ceeq=ceeq+1
ce,ceeq,(eps11*(a1*2))+(eps12*(a2*2))-(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1

nsel,s,node,,a1n_a2p_a3n_corner
*get,neg_node,node,0,num,min
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
nsel,a,node,,a1p_a2n_a3p_corner
pos_node=node(-x_,-y_,-z_)
ceeq=ceeq+1
ce,ceeq,(eps21*(a1*2))-(eps22*(a2*2))+(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1 	! x->2  
ceeq=ceeq+1
ce,ceeq,(eps31*(a1*2))-(eps32*(a2*2))+(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
ceeq=ceeq+1
ce,ceeq,(eps11*(a1*2))-(eps12*(a2*2))+(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1

nsel,s,node,,a1n_a2p_a3p_corner
*get,neg_node,node,0,num,min
	! 	get x,y,z locations of current node in active coord system
	x_=nx(neg_node)
	y_=ny(neg_node)
	z_=nz(neg_node)
nsel,a,node,,a1p_a2n_a3n_corner
pos_node=node(-x_,-y_,-z_)
ceeq=ceeq+1
ce,ceeq,(eps21*(a1*2))-(eps22*(a2*2))-(eps23*(a3*2)),neg_node,ux,-1,pos_node,ux,1 	! x->2  
ceeq=ceeq+1
ce,ceeq,(eps31*(a1*2))-(eps32*(a2*2))-(eps33*(a3*2)),neg_node,uy,-1,pos_node,uy,1	! y->3
ceeq=ceeq+1
ce,ceeq,(eps11*(a1*2))-(eps12*(a2*2))-(eps13*(a3*2)),neg_node,uz,-1,pos_node,uz,1	! z->1

allsel,all
*end