
DEBUT();

#read mesh

pmesh=LIRE_MAILLAGE(FORMAT='MED',);

pmodel=AFFE_MODELE(MAILLAGE=pmesh,
                   AFFE=_F(TOUT='OUI',
                           PHENOMENE='MECANIQUE',
                           MODELISATION='3D',),);

#define steel and apply to construction

kanagu=DEFI_MATERIAU(ELAS=_F(E=20600000.0,
                             NU=0.28,),);

hari=DEFI_MATERIAU(ELAS=_F(E=9e5,
                           NU=0.4,),);

matprops=AFFE_MATERIAU(MAILLAGE=pmesh,
                       AFFE=(_F(GROUP_MA='kanagu',
                                MATER=kanagu,),
                             _F(GROUP_MA='hari',
                                MATER=hari,),),);

#apply boundary conditions and load 

FLoad=AFFE_CHAR_MECA(MODELE=pmodel,
                     DDL_IMPO=(_F(GROUP_MA='kotei',
                                  DX=0.0,
                                  DY=0.0,
                                  DZ=0.0,),
                               _F(GROUP_MA='hanbun',
                                  DX=0,
                                  DY=0,),),
                     FORCE_FACE=_F(GROUP_MA='saika',
                                   FZ=-3.5e-1,),);

#                    PRES_REP =     (_F(GROUP_MA='topz',PRES=pressure,),),);
# calculate stresses with applied BC and loads:

result=MECA_STATIQUE(MODELE=pmodel,
                     CHAM_MATER=matprops,
                     EXCIT=_F(CHARGE=FLoad,),);
##==============postprocessing results - displacements etc



result=CALC_CHAMP(CRITERES=('SIEQ_ELNO',),
                  CONTRAINTE=('SIGM_ELNO','SIEF_ELGA',),
                  reuse =result,
                 MODELE=pmodel,
                 CHAM_MATER=matprops,
                 RESULTAT=result,
                 EXCIT=_F(CHARGE=FLoad,),);

#================

pmesh=DEFI_GROUP(reuse =pmesh,
                 MAILLAGE=pmesh,
                 CREA_GROUP_NO=_F(GROUP_MA='saika',
                                  NOM='saika',),
                 INFO=2,);

#                 CREA_GROUP_MA=_F(NOM='top',GROUP_MA=('top',),),
# calculate nodal equivalents and reaction forces


result=CALC_CHAMP(FORCE=('FORC_NODA',),
                    CRITERES=('SIEQ_NOEU',),
                    CONTRAINTE=('SIGM_NOEU',),
                    reuse =result,
               RESULTAT=result,);

# Create a table for displacements ('DEPL') and nodal forces ('FORC_NODA'). 
# In both cases the three #directions x, y and z (RESULTANTE=('DX','DY','DZ',)) will be extracted. 

TB_nodf=POST_RELEVE_T(ACTION=(_F(OPERATION='EXTRACTION',
                                 INTITULE='displacements',
                                 RESULTAT=result,
                                 NOM_CHAM='DEPL',
                                 TOUT_ORDRE='OUI',
                                 GROUP_NO='saika',
                                 RESULTANTE=('DX','DY','DZ',),
                                 MOYE_NOEUD='NON',),
                              _F(OPERATION='EXTRACTION',
                                 INTITULE='NodalForces',
                                 RESULTAT=result,
                                 NOM_CHAM='FORC_NODA',
                                 PRECISION=0.000001,
                                 GROUP_NO='saika',
                                 RESULTANTE=('DX','DY','DZ',),),),
                      TITRE='[Reaction] Nodal Forces',);

#write the data to a file (define a file in ASTK with file unit 26, in this case).

IMPR_TABLE(TABLE=TB_nodf,
           FORMAT='TABLEAU',
           UNITE=26,
           SEPARATEUR=' * ',
           TITRE='displacements at nodes on group TOP',);

# print results in med format (eg for Salome)

IMPR_RESU(FORMAT='MED',
          UNITE=21,
          RESU=_F(MAILLAGE=pmesh,
                  RESULTAT=result,
                  NOM_CHAM=('SIGM_NOEU','SIEQ_NOEU','SIEQ_ELNO','DEPL',),),);
#### now start geometrical non linear phase #####
#### -------------------------------------- #####


# determine stiffness matrix of the structure and DDL vector

Smech=CALC_MATR_ELEM(OPTION='RIGI_MECA',
                     MODELE=pmodel,
                     CHAM_MATER=matprops,
                     CHARGE=FLoad,);

# determine stress due to loading - create the geometrical stiffness matrix

Fsigma=CREA_CHAMP(TYPE_CHAM='ELGA_SIEF_R',
                  OPERATION='EXTR',
                  RESULTAT=result,
                  NOM_CHAM='SIEF_ELGA',
                  NUME_ORDRE=1,);

# linear model calculated                  
# calculate geometrical stiffness matrix due to loading of the structure

Sgeom=CALC_MATR_ELEM(OPTION='RIGI_GEOM',
                     MODELE=pmodel,
                     SIEF_ELGA=Fsigma,);

nddl=NUME_DDL(MATR_RIGI=Smech,);

SAmech=ASSE_MATRICE(MATR_ELEM=Smech,
                    NUME_DDL=nddl,);

SAgeom=ASSE_MATRICE(MATR_ELEM=Sgeom,
                    NUME_DDL=nddl,);

# multiply results by -1 or reverse force, since we solve 
# A*x = lamdba*B*x iso (A+mu*B)x=0

buckle=CALC_MODES(OPTION='PLUS_PETITE',
                        SOLVEUR_MODAL=_F(NMAX_ITER_SOREN=50,
                                         PREC_SOREN=2e-16,
                                         METHODE='SORENSEN',),
                        MATR_RIGI=SAmech,
                        MATR_RIGI_GEOM=SAgeom,
                        TYPE_RESU='MODE_FLAMB',
                        CALC_CHAR_CRIT=_F(NMAX_CHAR_CRIT=6,),
                        INFO=1,);

IMPR_RESU(FORMAT='MED',
          RESU=(_F(MAILLAGE=pmesh,
                   RESULTAT=buckle,
                   NOM_CHAM=('DEPL','SIGM_NOEU',),),
                _F(MAILLAGE=pmesh,
                   RESULTAT=result,
                   NOM_CHAM=('SIGM_ELNO','SIGM_NOEU','SIEQ_ELNO','SIEQ_NOEU','DEPL',),),),);

FIN();
