DEBUT();

mesh=LIRE_MAILLAGE(INFO=2,
                   UNITE=20,
                   FORMAT='MED',);

mesh=DEFI_GROUP(reuse=mesh,
                MAILLAGE=mesh,
                CREA_GROUP_NO=(_F(TOUT_GROUP_MA='OUI',),),);


model=AFFE_MODELE(MAILLAGE=mesh,
                  AFFE=_F(GROUP_MA=('stem',),
                         PHENOMENE='MECANIQUE',
                         MODELISATION='POU_D_T',),);

steel=DEFI_MATERIAU(ELAS=_F(E=210000.,NU=0.3,RHO=8.0e-9),);



material=AFFE_MATERIAU(
            MAILLAGE=mesh,
            AFFE=_F(GROUP_MA=('stem',),MATER=steel,),
			);
       
a=(1.0/8/10/2)**0.5*100;

elemcar=AFFE_CARA_ELEM(
   MODELE=model,
   POUTRE=_F(
         GROUP_MA=('stem',),
     SECTION='RECTANGLE',CARA=('HY','HZ',),VALE=(a,2*a),
	 ),
);


ground=AFFE_CHAR_MECA(MODELE=model,
                      DDL_IMPO=(
                          _F(GROUP_NO=('fix',),DX=0,DY=0,DZ=0,DRZ=0,),
                          _F(GROUP_NO=('endmass',),DX=0,DY=0,DRZ=0,),),);
                    
					#uncomment in a second run
					#to see the behavior with a fixed force
					#FORCE_NODALE=_F(GROUP_NO=('endmass',),FZ=-1000,),);

load=AFFE_CHAR_MECA(MODELE=model,
                    FORCE_NODALE=_F(GROUP_NO=('endmass',),FZ=-1000,),);
                      
#for variable load
resc11p1=MECA_STATIQUE(MODELE=model,
                       CHAM_MATER=material,
                       CARA_ELEM=elemcar,
                       EXCIT=(_F(CHARGE=ground,),
                              _F(CHARGE=load,),),
                       OPTION='SIEF_ELGA',);
#for fixed load
resc12p1=MECA_STATIQUE(MODELE=model,
                       CHAM_MATER=material,
                       CARA_ELEM=elemcar,
                       EXCIT=_F(CHARGE=ground,),
                       OPTION='SIEF_ELGA',);
#stress fields for variable load
sigc11p1=CREA_CHAMP(TYPE_CHAM='ELGA_SIEF_R',
                    OPERATION='EXTR',
                    RESULTAT=resc11p1,
                    NOM_CHAM='SIEF_ELGA',
                    TYPE_MAXI='MINI',
                    TYPE_RESU='VALE',);
#geometrical stiffness matrix for variable load
regc11p1=CALC_MATR_ELEM(OPTION='RIGI_GEOM',
                         MODELE=model,
                         CARA_ELEM=elemcar,
                         SIEF_ELGA=sigc11p1,);
#stress fields for fixed load
sigc12p1=CREA_CHAMP(TYPE_CHAM='ELGA_SIEF_R',
                    OPERATION='EXTR',
                    RESULTAT=resc12p1,
                    NOM_CHAM='SIEF_ELGA',
                    TYPE_MAXI='MINI',
                    TYPE_RESU='VALE',);
#geometrical stiffness matrix for fixed load
regc12p1=CALC_MATR_ELEM(OPTION='RIGI_GEOM',
                         MODELE=model,
                         CARA_ELEM=elemcar,
                         SIEF_ELGA=sigc12p1,);
#material stiffness matrix for both loads
remep1=CALC_MATR_ELEM(
      OPTION='RIGI_MECA',
      MODELE=model,
      CHAM_MATER=material,
      CARA_ELEM=elemcar,
      CHARGE=(ground,load,),
	  );
#matrix assemblies
nup1=NUME_DDL(MATR_RIGI=remep1,);

ramc1p1=ASSE_MATRICE(MATR_ELEM=remep1,NUME_DDL=nup1,);

ragep1=ASSE_MATRICE(MATR_ELEM=regc11p1,NUME_DDL=nup1,);

ragc12p1=ASSE_MATRICE(MATR_ELEM=regc12p1,NUME_DDL=nup1,);

#addition
ramep1=COMB_MATR_ASSE(
    COMB_R=(
	_F(MATR_ASSE=ramc1p1,COEF_R=1.0,),
    _F(MATR_ASSE=ragc12p1,COEF_R=1.0,),
	),
	);
#here we set the mini and maxi value of the modes with STURM
#the calculation fails if there are no modes in this range
#it is then necessary to adjust
#sturm u45201
#then we calucelate in the range
mini=-100;

maxi=100;

IMPR_STURM(
  MATR_A=ramep1,
  MATR_B=ragep1,
  TYPE_RESU='MODE_FLAMB',
  CHAR_CRIT_MIN=mini,
  CHAR_CRIT_MAX=maxi,
);

flamb=MODE_ITER_SIMULT(
    MATR_A=ramep1,
    MATR_B=ragep1,
    TYPE_RESU='MODE_FLAMB',
    CALC_FREQ=_F(
	   OPTION='BANDE',
       CHAR_CRIT=(mini,maxi,),
	   DIM_SOUS_ESPACE=80,
	   #NMAX_ITER_SOREN=80,
	   ),
);

flamb=NORM_MODE(reuse=flamb,
                MODE=flamb,
                NORME='TRAN',);


IMPR_RESU(FORMAT='MED',
          UNITE=80,
          RESU=_F(
		          RESULTAT=flamb,
                  NOM_CHAM='DEPL',
				  ),
);

IMPR_RESU(
    MODELE=model,FORMAT='RESULTAT',
    RESU=_F(RESULTAT=flamb,
                  TOUT_PARA='OUI',
                  FORM_TABL='OUI',
				  ),
);

FIN();
