/out
/com
/com *** LMATRIX Command Macro
/com
!/com Calculates incremental inductance matrix about an operating point.
!/com ANSYS Revision 5.5
!/com
!/com Theory based on the paper: M.Gyimesi, D.Ostergaard,
!/com "Inductance Computation by Incremental Finite Element Analysis",
!/com presented at the CEFC Conference, Tucson, AR, 1998.
!/com
!/com
!/com *************************************************
!/com
!/com ARGUMENTS
!/com  
!/com arg1 = Symmetry factor; defaults to 1
!/com arg2 = prefix for component names 
!/com arg3 = 1D array name for operating coil currents per turn
!/com arg4 = 2D array name for computed inductances, defaults to LMATRIX
!/com arg5 = help/print key
!/com        0 : print only computed inductance matrix - default
!/com        1 : print also current and energy increments
!/com        2 : print also fields due to incremented loads
!/com        3 : print everything, do not use scratch file
!/com
!/com USAGE
!/com
!/com Operating coil currents per turn must be described in 1D array ,
!/com defind in advance by *DIM,arg3,array,N , where N-number of coils. 
!/com Operating coil currents must be non-zero.
!/com Prescribe a negligible current for zero currents.
!/com Current density JS loads of field elements and
!/com current in real constant set of coil36 elements
!/com should reflect the number of turns.
!/com
!/com Coils must be incorporated in element components named coil1, coil2 ...
!/com unless a different prefix is provided by 2nd argument.
!/com The suffix of component names should match the pertinent coil.
!/com A coil component may incorporate different element types.
!/com
!/com Computed inductances are returned in 2D array: arg4(N,N).
!/com Computed inductances are printed in file, lmatrix.out
!/com
!/com RESTRICTIONS
!/com
!/com The macro works only with magnetic field elements: 53, 96, 97, 117.
!/com Further excitation may be described by current source element 36.
!/com Interface element, 115, can be applied to connect 96 and 97 regions.
!/com Infinite and circuit elements can also be used.
!/com However energy contribution from infinite elements is disregarded.
!/com
!/com Operating solutions must be obtained before calling LMATRIX
!/com by static analysis using frontal solver.
!/com
!/com Reserved file names:
!/com oper.esav, oper.emat, oper.tri, oper.db and lmatrix.out
!/com Make sure there is sufficient disk space available to store them.
!/com In case of a crash copy these files to file.* or jobname.*
!/com to recover operating solution files.
!/com
!/com The macro saves the database when called.
!/com If you wish to preserve an earlier database,
!/com save and rename it before calling the macro.
!/com
!/com The current in each coil36 component must be described by
!/com real contants which are not referred by other components.
!/com
!/com Standard ANSYS output file is used and switched at the end of macro.
!/com Macro exits at the BEGIN level and turns expanded print off (/nopr).
!/com
!/com Jobname can not be longer than 8 chararcters.
!/com
!/com *************************************************
!/com
!
/nopr                     ! turn off expanded printing
_ok3=0                    ! initialize for checking
_status=0
_esav=1                   ! files recording flags initialization
_emat=1
_tri=1
_db=1
_rmg=1
_rst=1
!
! ******** use scratch output unless debugging with arg5=3
!
*get,_cpu1,active,,time,cpu
_qq=arg5
*get,_ptyp,parm,_qq,type  ! obtain type of parameter
*if,_ptyp,ne,0,then
   *msg,error
   LMATRIX 5th argument must be numeric
   _ok3=1
*endif
!
!*if,arg5,le,2,then
!/out,scratch
!*endif
!
! ******** check jobname
!
*dim,_jname,char,2
*get,_jname(1),active,,jobnam,,start,1 ! get jobname
*get,_jname(2),active,,jobnam,,start,9 ! get jobname
!
!*stat,_jname
!
*if,_jname(2),ne,' ',then
*msg,error
LMATRIX is restricted to jobnames up to 8 characters.
_ok3=1
*endif
!
! ********* adjust symmetry factor
!
_qq=arg1
*get,_ptyp,parm,_qq,type  ! obtain type of parameter
*if,_ptyp,ne,0,then
   *msg,error
   Input value for symmetry factor must be numeric
   _ok3=1
*endif
!
*if,arg1,le,0,then
 _symm = 1.0
*else
 _symm = arg1
*endif
!
! ******** obtain prefix
!
_qq=
_qq=arg2
*get,_ptyp,parm,_qq,type  ! obtain type of parameter
*if,_ptyp,eq,0,then
   *msg,error
 Component name prefix must be alphanumeric string enclosed in single quotes
   _ok3=1
*else
   *if,_ptyp,eq,3,then
      _pre=_qq
   *else
      *msg,error
      Component name prefix must be an alphanumeric string
      _ok3=1
   *endif
*endif
!
! ******** obtain 1D array name
!
_qq=
_qq=arg3
*get,_ptyp,parm,_qq,type  ! obtain type of parameter
*if,_ptyp,eq,3,then
   *get,_ztyp,parm,arg3,type
   *if,_ztyp,ne,1,then
      *msg,error,_qq
       Coil current array %c must be defined before LMATRIX command
      _ok3=1
   *else
      _cur=arg3
   *endif
*elseif,_ptyp,eq,0,then
      *msg,error
 Coil current array name must be alphanumeric string enclosed in single quotes
      _ok3=1
*else
      *msg,error
        Coil current array name must be an alphanumeric string
      _ok3=1
*endif
!
! ******** obtain 2D array name
!
_qq=
_qq=arg4
*get,_ptyp,parm,_qq,type  ! obtain type of parameter
*if,_ptyp,eq,3,then
   *if,arg4,eq,' ',then
      arg4='lmatrix'
      *msg,note,arg4
      Inductance matrix array parameter name defaults to: %c
   *endif
*endif
*if,_ptyp,eq,0,then
   *if,_qq*1.0e+30,eq,0,then
      arg4='lmatrix'
      *msg,note,arg4
      Inductance matrix array parameter name defaults to: %c
   *else
      *msg,error
 Inductance matrix name must be alphnumeric string enclosed in single quotes
      _ok3=1
   *endif
*endif
!
cm,_elem,elem
*if,_ok3,ne,1,then
! ********* check analysis type
!
*get,_antyp,active,,anty
*if,_antyp,ge,1,then
    *msg,error
    LMATRIX macro is available only for static analysis
    _ok3=1
*endif
!
! ********* check declaration
! 
*get,_nrows,parm,_cur,dim,x   ! get number of rows, i.e. # of coils
*if,_status,gt,2,then
    *msg,error
    Non-existent or improperly dimensioned 1D-array
    _ok3=1
*endif

! **********  check valid elements
!
esel,u,ename,,53    ! valid element types
esel,u,ename,,96
esel,u,ename,,97
esel,u,ename,,117
esel,u,ename,,36
esel,u,ename,,115
*get,_mxvld,elem,,count
*if,_mxvld,gt,0,then
 *msg,error
   Some of the selected element types are not supported by LMATRIX
 _ok3=1
*endif
cmsel,s,_elem

*endif
!
! *********** declare working arrays and create component names
!
!                                ! declare working arrays
*dim,_ind,array,_nrows,_nrows     ! incremental inductance matrix
*dim,_ene,array,_nrows,_nrows    ! incremental energy matrix
*dim,_nami,char,_nrows           ! component names of coils
!
*do,_i,1,_nrows
_nami(_i)='%_pre%%chrval(_i)%'
*enddo ! of _i
!
!******************* get further active data
!
*get,_mnu3,active,,menu       ! current menu
*get,_level,active,,rout      ! current routine level
*get,_neqit,active,,solu,eqit ! current number of equilibrium iterations
!
/out,scratch
finish                        ! go to begin level
/out
!
*if,_ok3,ne,1,then
/gopr
/out,data
! *********** check if components have elements
!
_noelm=0
!
*do,_i,1,_nrows
cmsel,s,_nami(_i)
*stat,_nami  !  12/26
*if,_status,gt,2,then
 *msg,error
 Coil element components are not properly defined; check 2nd argument
 _ok3=1
 *exit
*endif
*get,_ecount,elem,,count! get number of elements
*if,_ecount,lt,1,then
 _noelm=1
*endif
*enddo ! of _i
!
*if,_noelm,eq,1,then
 *msg,error
 Coil component contains no elements
 _ok3=1
*endif

*endif
!
cmsel,s,_elem
!
! ************ save operating solutions ***************
!
*if,_ok3,ne,1,then
 save
 /rename,%_jname(1)%,esav,,oper,esav
 _esav=
 *if,_status,gt,2,then
  _esav=1
  *msg,error,_jname(1)
  File %c.esav is required for LMATRIX command
  _ok3=1
 *endif
 /rename,%_jname(1)%,emat,,oper,emat
 _emat=
 *if,_status,gt,2,then
  _emat=1
  *msg,error,_jname(1)
  File %c.emat is required for LMATRIX command
  _ok3=1
 *endif
 /rename,%_jname(1)%,tri,,oper,tri
 _tri=
 *if,_status,gt,2,then
  _tri=1
  *msg,error,_jname(1)
  File %c.tri is required for LMATRIX command
  _ok3=1
 *endif
 /rename,%_jname(1)%,db,,oper,db
 _db=
 *if,_status,gt,2,then
  _db=1
  *msg,error,_jname(1)
  File %c.db is required for LMATRIX command
  _ok3=1
 *endif
/out,scratch
/uis,msgpop,3
/nerr,,,-1
/uis,msgpop
/out
/copy,%_jname(1)%,rmg,,oper,rmg
  _rmg=
  _rst=1
 *if,_status,gt,2,then
  _rmg=1
  _rst=0
 /copy,%_jname(1)%,rst,,oper,rst
  *if,_status,gt,2,then
  _rst=1
  /nerr
   *msg,error
   Results file for Operating point solution is missing
   _ok3=1
  *endif
 *endif
*endif
!
! ************ compute nominal and incremental current densities
!

*if,_ok3,ne,1,then

_curmax=0.0
!
 *do,_i,1,_nrows
 *if,abs(_cur(_i)),gt,_curmax,then
 _curmax=abs(_cur(_i))
*endif
*enddo
!
*if,_curmax,lt,1.0e-6,then
 *msg,error
   Nominal applied currents from operating point solution are zero
  _ok3=1
*endif
!
_tinycur=1.0e-6*_curmax
!
*do,_i,1,_nrows
!
*if,abs(_cur(_i)),lt,_tinycur,then
 *msg,error
  Current in one of the coils is zero, nonzero current required
 _ok3=1
*endif
*enddo

*endif
!
! ************** end of checking
! 
*if,_ok3,ne,1,then
! 
! ************* field analyses with increment currents in coils, _ii and, _jj
!
*do,_ii,1,_nrows
*do,_jj,_ii,_nrows
parsav,all                       ! save parameters
!
!/com
!/com linearized incremental solutions with no save
!/com
!!!                              ! restore operating solution
/copy,oper,esav,,%_jname(1)%,esav
/copy,oper,emat,,%_jname(1)%,emat
/copy,oper,tri,,%_jname(1)%,tri
/copy,oper,db,,%_jname(1)%,db       ! resume from operating DB
resume
parres                           ! restore parameters
!
/solu
!
! ------------------ set up linearized incremental analysis
!
!                           ! linearized incremental analysis
magopt,4,,1                 ! Incremental solution with no save
antype,static,rest          ! static restart analysis from operating point
kuse,1                      ! linearized
neqit,1                     ! force 1 eq iteration otherwise NR diverges
eqslv,front                 ! apply front solver
autots,off                  ! turn off auto time stepping
nlgeom,off                  ! turn off large deformations
nsubs,1                     ! number of substeps is one
!
! ------------------  scale according to incremental loads
!
! --- scale down everything
!
_scale=1.0e-12
bfescal,JS,_scale
!
! --------------- scale real of 36 elements
!
*get,_nele,elem,,count             ! number of selected elements
!*stat,_nele
!
*do,_iele,1,_nele                  ! loop on selected elements
!*stat,_iele
!
*get,_e1,elem,,num,max             ! get largest elem # remaining in the set
!*stat,_e1
!
*if,_e1,eq,0,exit
!
*get,_ityp,elem,_e1,attr,type      ! obtain type number
*get,_rcon,elem,_e1,attr,real      ! obtain real number
*get,_jtyp,etyp,_ityp,attr,enam    ! obtain stiffness number of _ityp
!
!*stat,_ityp
!*stat,_jtyp
!*stat,_rcon
!
*if,_jtyp,eq,36,then
!
*get,_rcur,rcon,_rcon,const,2      ! 2nd real is current for coil36
!*stat,_rcur
rmodif,_rcon,2,_scale*_rcur        ! scale current in real set for coil36
esel,u,real,,_rcon                 ! unselect 36 types with processed real
!
*else
esel,u,type,,_ityp                 ! unselect not-36 types in the component
*endif
!
!elist
!
*enddo  ! of _iele; --- end of real scale 36 elements
!
! ---- end of scale down
!
!/com ### after down scale
!bfelist
!rlist
!
*do,_kkk,1,_nrows          ! loop on coils
!
cmsel,s,_nami(_kkk)
!*stat,_kkk
!elist
!
*if,_kkk,eq,_ii,then
*if,_kkk,eq,_jj,then
!                           ! kkk=ii and kkk=jj
_scale=1.0e+12
*else
!                           ! kkk=ii and kkk<>jj
_scale=1.0e+12
*endif
*else
*if,_kkk,eq,_jj,then
!                           ! kkk<>ii and kkk=jj
_scale=1.0e+12
*else
!                           ! kkk<>ii and kkk<>jj
_scale=1.0
*endif
*endif
!
_scale=_scale*_curmax/_cur(_kkk)
!
bfescal,JS,_scale                  ! scale JS body force of input DB
!
! --------------- scale real of 36 elements
!
*get,_nele,elem,,count             ! number of selected elements
!*stat,_nele
!
*do,_iele,1,_nele                  ! loop on selected elements
!*stat,_iele
!
*get,_e1,elem,,num,max             ! get largest elem # remaining in the set
!*stat,_e1
!
*if,_e1,eq,0,exit
!
*get,_ityp,elem,_e1,attr,type      ! obtain type number
*get,_rcon,elem,_e1,attr,real      ! obtain real number
*get,_jtyp,etyp,_ityp,attr,enam    ! obtain stiffness number of _ityp
!
!*stat,_ityp
!*stat,_jtyp
!*stat,_rcon
!
*if,_jtyp,eq,36,then
!
*get,_rcur,rcon,_rcon,const,2      ! 2nd real is current for coil36
!*stat,_rcur
rmodif,_rcon,2,_scale*_rcur        ! scale current in real set for coil36
esel,u,real,,_rcon                 ! unselect 36 types with processed real
!
*else
esel,u,type,,_ityp                 ! unselect not-36 types in the component
*endif
!
!elist
!
*enddo  ! of _iele; --- end of real scale 36 elements
!
*enddo ! of coils, kkk
!
! --------------- end of scaling; solve and obtain energies
!
cmsel,s,_elem
!
!/com ### after scale
!bfelist
!rlist
!
biot,new
!
outres,esol,all
/uis,msgpop,4
solve
/uis,msgpop
*if,_status,gt,2,then,exit
     *msg,error
  LMATRIX incremental solution encountered an error condition
  _ok3=1
*endif
!-------------------- restore settings 
/out,scratch
finish
/out
/solu
outres
neqit,_neqit
kuse
/out,scratch
fini
/out
!
/post1
!              
!etab,_ene,SMISC,5
etab,_ene,IENE
/out,scratch
ssum
/out
*get,_ene(_ii,_jj),SSUM,,ITEM,_ene
_ene(_jj,_ii)=_ene(_ii,_jj)
_eneiijj=_ene(_ii,_jj)
!
*if,arg5,ge,1,then         ! print increment details
/out
/com
/com linearized incremental solution from operating point
/com
*if,_ii,eq,_jj,then
*vwrite,_jj,_curmax
('current incremented in coil',f3.0,' by',e15.6)
*else
*vwrite,_ii,_jj,_curmax
('current incremented in coils',f3.0,' and',f3.0,' by',e15.6)
*endif
!
*if,arg5,ge,2,then         ! print even field details
etab,_x,cent,x
etab,_y,cent,y
etab,_z,cent,z
etab,_hsum,h,sum
etab,_bsum,b,sum
pretab,_x,_y,_z,_hsum,_bsum,_ene
!presol,h
!presol,b
!
*endif
!
*vwrite,_eneiijj
(' total energy increment=',e15.6)
!
!*if,arg5,le,2,then
!/out,scratch
!*endif
!
*endif ! of print incremental details: arg5 >= 1 
/out,scratch
fini
/out
!
*enddo ! of jj coil excitation increment
*enddo ! of ii coil excitation increment
!
! ************* end of field analyses with incremented currents
!
! ************* compute L matrix
!
*do,_ii,1,_nrows
*do,_jj,_ii,_nrows
*if,_ii,eq,_jj,then
!                     ! self inductance
_ind(_ii,_jj)=_ene(_ii,_ii)*2.0/_curmax**2*_symm
*else
!                     ! mutual inductance
_ind(_ii,_jj)=(_ene(_ii,_jj)-_ene(_ii,_ii)-_ene(_jj,_jj))/_curmax**2*_symm
_ind(_jj,_ii)=_ind(_ii,_jj)
*endif
!
*enddo ! of jj coil excitation increment
*enddo ! of ii coil excitation increment
!
! ************ generate output
!
*get,_cpu2,active,,time,cpu
 _dcpu=_cpu2-_cpu1
/nopr
/out,lmatrix,out
 *msg,info
%/____________ INDUCTANCE CALCULATION SUMMARY ___________________

/com
!
 *do,_ii,1,_nrows
 _indii=_ind(_ii,_ii)
 *vwrite,_ii,_indii
 ('Self inductance of coil ',f3.0,' = ',e15.5,' (H)')
 *enddo
!
 *do,_ii,1,_nrows
 *do,_jj,_ii+1,_nrows
 _indij=_ind(_ii,_jj)
 *vwrite,_ii,_jj,_indij
 ('Mutual inductance between coils ',f3.0,' and ',f3.0,' = ',e15.5,' (H)')
 *enddo
 *enddo
!
 *vwrite,arg4,_nrows,_nrows
 (/'Inductance matrix is stored in array parameter ',a,'(',f3.0,',',f3.0,')')
 /com
 *vwrite,arg4
 ('Inductance matrix is stored in file ',a,'.txt'/)
 *msg,info
%/_____________________________________________________________
*vwrite,_dcpu
 ('Time (cpu):     ',f6.1,'Sec')

 /com
!
/out
/gopr
!
! -------------- define output inductance array with name 'arg4'
!
*set,%arg4%      !   3/99
*dim,%arg4%,array,_nrows,_nrows
*do,_ii,1,_nrows
*do,_jj,1,_nrows
  %arg4%(_ii,_jj,1)=_ind(_ii,_jj,1)
*enddo
*enddo
!
 *if,_mnu3,ne,1,then     ! menu off
  *list,lmatrix,out      ! batch listing
  *else
  *list,lmatrix,out      ! batch listing
  *if,_mg1,eq,0,then
    *uilist,lmatrix,out    ! interactive listing
  *endif
 *endif
!
! -------------- write inductance matrix array in output file
!
/nopr
/out,%arg4%,txt
/com
/com ****** COMPUTED INDUCTANCE ARRAY ******
/com
/com        I       J         Value
/com
!
*do,_ii,1,_nrows
*do,_jj,_ii,_nrows
 _indij=_ind(_ii,_jj)
 *vwrite,_ii,_jj,_indij
 (5X,2(F5.0,3X),E13.5)
*enddo
*enddo
!
/com
/com ****** END OF INDUCTANCE ARRAY ****** 
/com
!
/out
!
! -------------- reinstall operating database
!
!

*endif  ! of _ok3
parsav,all                       ! save parameters
/out,scratch
/solu
magopt,4                         ! clear magnsv lopgen/soptcm/edgclr
fini
*if,_esav,ne,1,then
/rename,oper,esav,,%_jname(1)%,esav
*endif
*if,_emat,ne,1,then
/rename,oper,emat,,%_jname(1)%,emat
*endif
*if,_tri,ne,1,then
/rename,oper,tri,,%_jname(1)%,tri
*endif
*if,_db,ne,1,then
/rename,oper,db,,%_jname(1)%,db
*endif
*if,_rmg,ne,1,then
/rename,oper,rmg,,%_jname(1)%,rmg
*endif
*if,_rst,ne,1,then
/rename,oper,rst,,%_jname(1)%,rst
*endif

/out

*if,_ok3,ne,1,then
resume
*endif
parres                           ! restore parameters
!
! 
! ************************** clean up
!
cmdele,_elem
_ind =
_nrows= $_mnu3= $_ok3= $_symm= $_cchk3= $_jname=
_status= $_ene= $_curmax= $_tinycur= $_absi= $_scale=
_indii= $_indij= $_nami= $_ecount= $_cur= $_noelm= $_ztyp=
_antyp= $_level= $_neqit= $_ityp= $_jtyp= $_rcon= $_rcur= $_e1=
_ptyp= $_qq= $_pre= $_mxvld= $_ind= $_eneiijj= $_nami=
_nele= $_iele= $_esav= $_emat= $_tri= $_rmg= $_rst= $_db=
_x= $_y= $_z= $_hsum= $_bsum= $_ii= $_jj= $_kkk= $_i=
!
! menu, neqit, kuse, rutine, magopt, real, js, nopr
!
/out