$ FILE MMECHA.V707 $ MASSLESS MECHANISM DETECTOR FOR REAL MODES
$ FIND MASSLESS MECHANISM MODE SHAPES.  MAG  OCTOBER 28, 1999
$ CONSTRAIN MASSLESS MECHANISM MODES AUTOMATICALLY
$ experimental version tuned for Ford, mag, March 6, 2000
$ put mtiny on matgpr statements
$ updated 24 Mar 00.  LaGrange multipliers on mechansisms replaced
$ by spc-type operations.  Mechanisms with mass not constrained.  Changes in 
$ lower case.  Mag
$ updated April 5, 2000.  Generalized damping of mechanism modes added, GMA.  
$ Changes in lower case
$ECHOOFF $

COMPILE SEMG NOLIST NOREF $ 
ALTER 123 $ AFTER EMA
CALL DBSTORE GPLS,,,,//0/0/'     '/0 $ GPLS NEEDED FOR MATGPR

COMPILE SEMODES NOLIST NOREF $
ALTER 'LL MODEFS' $
CALL LOADCHK USET,GOT,GM,,GPLS,SILS/SCRQ $

COMPILE SEMFREQ NOLIST NOREF $ SOL 111
$ COMPUTE THE GENERALIZED FORCE OF THE MECHANISM MODES
ALTER 'FRLG ' $ AFTER DYNAMIC LOADS CALCULATEDT
CALL LOADCHK USETD,GOD,GMD,PPF,GPLD,SILD/QDG $ QDG, LOADS FROM AUTOCONST.
ALTER 'CALL FREQRS' $
CALL XADD5 QDG,QSGF,,,/QSGX/CU/CU/CZ/CZ/CZ $
EQUIVX QSGX/QSGF/ALWAYS $

COMPILE SEMTRAN NOLIST NOREF $
$ COMPUTE THE GENERALIZED FORCE OF THE MECHANISM MODES
ALTER 'TRLG ' $ AFTER DYNAMIC LOADS CALCULATED
CALL LOADCHK USETD,GOD,GMD,PPT,GPLD,SILD/QXG
CALL DBSTORE QXG,,,,//0/0/'     '/0 $

Compile gma nolist noref $
Alter 334 $ fluids, k4 not attempted
TYPE PARM,,RS,Y,MECHFIL=1.E-10 $ FILTER FOR MECHANISM MATRIX
Type db gpls,sils $ fetch by type statement
If (nob2ss>0 and not ue) then $
  Call dbfetch /mecha,,,,/0/0/0/0/0 $
  Smpyad mecha,b2ss,mecha,,,/genbb/3////1 $
  Diagonal genbb/genb/'column' $
  Norm genb/scrb///s,n,maxgenb $
  If (maxgenb>mechfil) then $
    Message //'uwm.  Generalized damping diagonal'/
    ' terms for deleted mechanisms.' $
    message //'These dampers not loaded for modal solutions' $
    Matprn genb// $
    Mpyad b2ss,mecha,/fdamp $
    MATGPR GPLs,USETD,SILs,fdamp//'H'/'d'//mechfil $
  Endif $
Endif $

COMPILE DISPRS NOLIST NOREF $
ALTER 85 $ 
TYPE PARM,,CS,N,CZ,CU=(1.,0.) $
IF (APP='TRANRESP') THEN $
  CALL DBFETCH /QXG,,,,/0/0/0/0/0 $
  TYPE PARM,,I,N,NT,NTM1 $
  PARAML QXG//'TRAILER'/1/S,N,NT $
  IF (NT>0) THEN $
    NTM1 = NT - 1 $
    MATGEN ,/VTM1/6/NT/NTM1/1 $
    PARTN QXG,VTM1,/QXGM1,,,/0 $
    CALL XADD5 QPF,QXGM1,,,/QPFX/CU/CU/CZ/CZ/CZ $
    EQUIVX QPFX/QPF/ALWAYS $
  ENDIF $
ENDIF $

COMPILE LOADCHK NOLIST NOREF $
SUBDMAP LOADCHK USETD,GOD,GMD,PPF,GPLD,SILD/QXG $
$ PRINTS LOADS THAT WILL BE DISCARDED BECAUSE OF MASSLESS MECH.
$ PUTS DISCARDED LOADS IN QXG FOR SPCFORCE PRINT
TYPE PARM,,I,N,NM,NF $ 
TYPE PARM,,CHAR8,N,OUTSET='P    ' $
CALL DBFETCH /MECHA,,,,/0/0/0/0/0 $
PARAML MECHA//'TRAILER'/1/S,N,NM $ NUMBER OF MECHANISMS
IF (NM=0) RETURN $
PARAML PPF //'TRAILER'/1/S,N,NF $ COLUMNS OF P
IF (NF=0) OUTSET = 'G    ' $
IF (OUTSET='G    ') RETURN $
SDR1  USETD,,MECHA,,,GOD,GMD,,,,/MECHPM,,/1/'STATICS' $
IF (NM>1) THEN 
  MATGEN, /UM/6/NM/0/NM $
  MPYAD MECHPM,UM,/V01P $ VECTOR THAT IS UNION OF ALL POINTS ON
  $ ALL MECHANISMS
ELSE $
  EQUIVX MECHPM/V01P/ALWAYS $
ENDIF $
$MESSAGE //'UIM.  UNION OF ALL POINTS IN ALL MECHANISMS' $
$MESSAGE //'ANY LOADS APPLIED TO THESE POINTS GO TO GROUND, NOT'/
$  ' TO THE STRUCTURE.' $
$ allow suppression of part or all of output.  Set mtiny large to
$ print little or no data
$type parm,,rs,y,mtiny=0.9 $ 
$message //'param, mtiny has a value of '/mtiny/
$  '. set it to a smaller value to see more DOFs listed' $
$MATGPR GPLD,USETD,SILD,V01P//'H'/OUTSET//mtiny $
$IF (OUTSET='G    ') RETURN $
PARTN PPF,,V01P/,QX1,,/0 $
MERGE, ,QX1,,,,V01P/QXG1/0 $
ADD QXG1,/QXG/-1.0 $
MPYAD MECHPM,PPF,/GENMF/1 $ M ROWS, F COLUMNS
MATGEN ,/UF/6/NF/0/NF $
MPYAD GENMF,UF,/GENMAX $
NORM GENMF/SCRN///S,N,MAXGENF
MESSAGE //'MAXIMUM GENERALIZED FORCE IN MECHANISM MODES'/MAXGENF $
IF (MAXGENF>1.E-10) THEN $
  MESSAGE //'USER WARNING ERROR, MASSLESS MECHANISMS ARE LOADED' $
  MESSAGE //'LARGEST GENERALIZED FORCE FACTOR IS '/MAXGENF $
  PARAML GENMAX//'TRAILER'/2/S,N,NN $
  IF (NN>1) THEN $
    MESSAGE //'GEN. FORCE FOR EACH MASSLESS MECHANISM IS' $
    MATPRN GENMAX// $
  ENDIF $
ENDIF $
RETURN $
END $ LOADCHK

COMPILE MODERS NOLIST NOREF $ CALLED FROM SOLS 103, 111, 112, 14X
ALTER 'READ.*KXX','' $
TYPE PARM,,I,N,SID=-2 $ SAYS USE THE EIGRL FIELD, BUT UPDATE IT WITH METHOD, 
CALL READX KXX,MXX,MR,DXR,EED,USET,CASES,VAXF,SILS,,
         ,LLL,EQEXINS/         
         LAMA,PHIX,MIX,OEIGS,EIGVMAT,,MECHA/                              
         'MODE '/S,NEIGV/NSKIP/SECND/'A    '/SID/'LAN  ' $ METHC DISABLED                              
CALL DBSTORE MECHA,,,,//0/0/'     '/0 $

COMPILE READX NOLIST NOREF $
SUBDMAP READX     KXX,MXX,MR,DXR,DYNAMICS,USET,CASES,
         VAXF,SILS,USET,DUM111,LLL,EQEXINS/         
         LAMA,PHIX,MIX,CMOEIGS,ROOTS,DUM214,MECHA/                              
          READAPP/NEIGV/NSKIP/METHF/CMSET/SID/METHOD $ METHC DISABLED                              
$ DUMI ARE NOT READ IN MODERS, V70.7
TYPE PARM,,LOGICAL,N,MMECH,METHF $
TYPE PARM,,CHAR8,N,READAPP,METHOD,CMSET $
TYPE PARM,,I,N,NEIGV,NSKIP$
TYPE PARM,,CS,N,CU=(1.0,0.0),CZ $ UTILITY PARAMS
TYPE PARM,,I,N,SID $ =-2 SAYS USE THE EIGRL FIELD, BUT UPDATE IT WITH METHOD, 
CALL MMFIND KXX,MXX,USET,SILS,EQEXINS,VAXF/
   MECHA,MECHLD,R/S,MMECH $
IF (MMECH) THEN $
  MESSAGE //'UWM. ASET DOFS IN MECHA ACTIVE IN MASSLESS MECHANISMS.' $
  MESSAGE //'EACH COLUMN IS A SEPARATE MECHANISM.  '/
            'LOADS ON THESE DOFS ARE DISCARDED'  $
  CALL DBFETCH /GPLS,,,,/0/0/0/0/0 $
type parm,,rs,y,mtiny=0.9 $ 
message //'param, mtiny has a value of '/mtiny/
  '. set it to a smaller value to see more DOFs listed' $
  MATGPR GPLS,USET,SILS,MECHA//'H'/'A'//mtiny $ 
  TYPE PARM,,CS,Y,EPSP=(1.E-10,0.0) $  
  TYPE PARM,,CHAR8,Y,DEBUG='NO' $
  IF (DEBUG='YES') THEN 
    MATPRN R// $
    DIAGON(30) $
    MATPRN KXX// $
    DIAGOFF(30) $
  ENDIF $
  SID = -2 $ PUT IN A LONGER READ LIST LATER
$ replace lmt with spcs
$  PARAML MECHLD//'TRAILER'/1/S,N,NMECHS $
$  PARAML KXX//'TRAILER'/1/S,N,NX $
$  TYPE PARM,,I,N,NC $
$  NC = NX + NMECHS
$  MATGEN, /VCXM/6/NC/NX/NMECHS $
$  TRNSP MECHLD/MECHLDT $
$  MERGE KXX,MECHLDT,MECHLD,,VCXM,/KXX2/-1 $
$  MERGE MXX,,,,VCXM,/MXX2/-1 $
  Paraml mechld//'trailer'/1/s,n,nmechs $ number of mechanisms
  If (nmechs>1) then $
    Matgen ,/un/6/nmechs/0/nmechs $
    Mpyad mechld,un,/vspc $ put all terms in one column.
  Else $
    Equivx mechld/vspc/always $
  Endif $ 1.0 where mechanisms to be grounded. 
  partn kxx,vspc,/kxx2,ksx,kxs,kss/-1 $
  partn mxx,vspc,/mxx2,msx,mxs,mss/-1 $
ELSE $
  EQUIVX KXX/KXX2/ALWAYS $
  EQUIVX MXX/MXX2/ALWAYS $
ENDIF $
READ     KXX2,MXX2,MR,DXR,DYNAMICS,USET,CASES,VAXF,SILS,,,LLL,EQEXINS/  
         LAMA,PHIX,MIX,OEIGS,ROOTS,/ 
         READAPP/S,N,NEIGV/NSKIP/METHF/CMSET/SID/METHOD $
IF (MMECH) THEN $
$  PARTN PHIX,,VCXM/PHIXX,QX,,/0 $ DISCARD QX, COMP. ZEROS
$  EQUIVX PHIXX/PHIX/ALWAYS $ 
  merge phix,,,,,vspc/phixx/0 $ insert nulls for spcs
  equivx phixx/phix/always $
  $ qx not being generated for now
ENDIF $
RETURN $
END $ READX

COMPILE MMFIND NOLIST NOREF $
SUBDMAP MMFIND K,M,USET,SILS,EQEXINS,VAWX/
  MECHU,MECHLD,R/MMECH $
$ MECHU IS THE MECHANISM MODES, ONE COLUMN PER MECHANISM
$ MECHLD HAS ONE NON-ZERO TERM PER COLUMN, USED AS CONSTRAINTS LATER
$ IF THERE ARE NONE MECHU IS PURGED AND MMECH=FALSE
$ R IS A VECTOR WITH UNIT TERMS AT THE END OF MECHANISMS.
$ MECHU IS NM VECTORS, ONE FOR EACH MECHANISM MODE, NORMALIZED
TYPE PARM,,LOGICAL,N,MMECH $ TRUE WHEN MASSLESS MECHANISMS EXIST
TYPE PARM,,RS,Y,SHIFT1=-1.234, MAXRATX $ USER MAY CHANGE TEST VALUE
TYPE PARM,,CS,N,CSHIFT1 $ COMPLEX FORM FOR ADD
TYPE PARM,,CS,N,CU=(1.0,0.0),CZ $ UTILITY PARAMS
TYPE PARM,NDDL,RS,Y,MAXRATIO $ NDDL DEFAULT FOR V71 = 1.E7
MAXRATX = MAXRATIO $ LOCAL VALUE
IF (MAXRATIO>1.E10) THEN $
  MAXRATx = 1.E10 $
  MESSAGE //'USER INFORMATATION MESSAGE.  MAXRATIO SET TO AN EXCSSIVE'/
   ' VALUE.  IT HAS BEEN RESET' $
ENDIF $
MESSAGE //'UIM.  MAXRATIO VALUE USED for singularity control IS'/
  MAXRATX $ 
PARAML M//'PRESENCE'////S,N,NOM $ -1 WHEN PURGED
CSHIFT1 = -CMPLX(SHIFT1) $ SHIFT IS DEFINED AS K-SHIFT1*M
CALL XADD5 K,M,,,/KM1/CU/CSHIFT1/CZ/CZ/CZ $
DCMP USET,SILS,EQEXINS,KM1,VAWX,/LM1,,////MAXRATX/'A'//26 $
MATMOD LM1,,,,,/T,LP1/21 $ TRID. TERMS, LOWER TRI. (L), PERMUTATION(P)
$ THE OPTION ABOVE WORKS ONLY AFTER V70.5
TYPE PARM,,LOGICAL,N,RAT1 $
CALL HIGHRAT KM1,T,LP1/ABSRAT1/S,RAT1 $ 1.0 AT END OF MASSLESS MECH.
NORM ABSRAT1/SCR1///S,N,MAXRAT1 $
$DIAGONAL ABSRAT1/R/'WHOLE'/0.0 $ SET NON-ZERO TERMS TO UNITY.  WATCH
$ OUT FOR PACKED ZEROS.
IF (MAXRAT1<MAXRATX) THEN $
  RAT1 = FALSE $
ENDIF $
IF (RAT1) THEN $
  EQUIVX ABSRAT1/MECHLOC/ALWAYS $ NON-ZERO AT END OF MECHANISMS
  MATMOD MECHLOC,,,,,/MECHLOCC,/28 $ DIAGONALIZE
  PARTN MECHLOCC,MECHLOC,/,,MECHLDX,/0 $ DELETE NULL COLUMNS
  NORM MECHLDX/MECHLD $ NORMALIZE TO AVOID OVERFLOWS
  MMECH = TRUE
  FBS LM1,,MECHLD/MECHU1 $ EACH COLUMN HAS ALL POINTS IN ONE MECHANISM
  NORM MECHU1/MECHU2 $ MAKE LARGEST TERMS 1.0
  Mpyad m,mechu2,/inertia $ mass loads in potential massless mechs
  Add mechu2,inertia/xmx///1 $ element multiply
  Paraml xmx//'trailer'/2/s,n,msize $
  Matgen ,/umsize/6/msize/0/msize $ unit vector
  Mpyad xmx,umsize,/genmass/1 $ diagonal terms of gen. Mass for mech modes
  TYPE PARM,,RS,Y,MECHFIL=1.E-10 $ FILTER FOR MECHANISM MATRIX
  Matmod genmass,,,,,/vmmechs,/2////mechfil $ discard numerical noise
  Paraml vmmechs//'trailer'/5/s,n,nzwdsm
  Paraml vmmechs//'trailer'/4/s,n,mtype $
  Paraml vmmechs//'trailer'/2/s,n,rowsmm
  Nzwdsm = nzwdsm/mtype $
  If (nzwdsm<rowsmm) then $
    Partn mechu2,vmmechs,/mechu3,,,/0 $
    Partn mechld,vmmechs,/mechld1,,,/0 $
    Equivx mechld1/mechld/always $
  Else $
    Equivx mechu2/mechu3/always $
  Endif $
  MATMOD MECHU3,,,,,/MECHU,/2////MECHFIL $ DISCARD NUMERICAL NOISE
ENDIF $
Paraml mechu//'presence'////s,n,nomechu $
If (nomechu=-1) mmech = false $
RETURN $
END $ mmfind

COMPILE HIGHRAT NOLIST NOREF $
SUBDMAP HIGHRAT KM,T,LP/ABSRAT/RAT $ 1.0 AT END OF MASSLESS MECH.
$ ABSRAT= FILTERED ABS(T/DIAG(KM)), KM=K+s*M
TYPE PARM,,LOGICAL,N,RAT $ TRUE IF ABSRAT NOT NULL. 
TYPE PARM,,I,N,NT2 $
TYPE PARM,NDDL,RS,Y,MAXRATIO $ NDDL DEFAULT FOR V71 = 1.E7
TYPE PARM,,RS,N,MAXRATX $ LOCAL VALUE
MAXRATX = MAXRATIO $
IF (MAXRATIO>1.E10) MAXRATx = 1.E10 $
DIAGONAL KM/KMDD/'SQUARE' $ DIAGONAL TERMS ALONG DIAGONAL
PARAML T//'TRAILER'/1/S,N,NT $ SIZE OF PROBLEM
NT2 = NT*2 $
MATGEN ,/VT2/6/NT2/NT/NT $ USED TO TAKE APART PACKED L, P MATRICES
PARTN LP,VT2,/,,P,/0 $ PERMUTATION MATRIX
SMPYAD P,T,P,,,/T1/3////1 $ RESEQUENCE TO EXTERNAL
DIAGONAL T1/TDD/'SQUARE' $ ALL BUT TWO BY TWOS REPRSENTED HERE
ADD T1,TDD/TOFFDD//-1.0 $ TWO BY TWOS HERE
PARAML TOFFDD//'TRAILER'/5/S,N,NZTOFF $ NON-ZERO TERMS IN OFF-DIAGS
IF (NZTOFF>0) THEN $
  MATGEN ,/UT/6/NT/0/NT $ UNIT VECTOR
  MPYAD TOFFDD,UT,/TOFFD $ VECTOR TO BLOCK 2 BY 2 PROCESSING
  MATMOD TOFFD,,,,,/TOFFDD1,/28 $ DIAGONALIZE
  ADD TOFFDD1,KMDD/KMDD1///3 $ 1.0 ON 2 BY 2 DIAGONAL TERMS THAT ARE NULL
  NORM TDD/SCRT///S,N,MAXTDD $ THIS BLOCKS THEM BEING IN HIGH MAXRATIOS
  IF (MAXTDD<1.E6) MAXTDD=1.E6
  TYPE PARM,,CS,N,CMAXTDD $
  CMAXTDD = CMPLX(MAXTDD) $
  ADD TOFFDD1,TDD/TDD1/CMAXTDD $ ADD LARGE TERMS TO DIAGONALS OF 2 BY 2'S
ELSE $
  EQUIVX KMDD/KMDD1/ALWAYS $ 
  EQUIVX TDD/TDD1/ALWAYS $
ENDIF $
DIAGONAL KMDD1/KMD/'COLUMN' $
DIAGONAL TDD1/TD/'COLUMN' $
ADD KMD,TD/RAT///2 $ ELEMENT DIVIDE
DIAGONAL RAT/ABSRAT1/'WHOLE'/1.0 $ MAKE ALL SIGNS PLUS
NORM ABSRAT1/SCR1///S,N,MAXRAT $
RAT = FALSE
IF (MAXRAT<MAXRATX) RETURN $
MATMOD ABSRAT1,,,,,/ABSRAT,/2////MAXRATX $ KEEP ONLY LARGE TERMS
TYPE PARM,,I,N,NZRAT $ GUARD AGAINST PURGED ABSRAT
PARAML ABSRAT//'TRAILER'/5/S,N,NZRAT $
IF (NZRAT>0) RAT = TRUE $
RETURN $
END $ MMECHA ALTER PACKAGE
ECHOON $