* fichier :  nafems-le3.dgibi

*** Options ...

    opti dime 3 mode trid elem seg2 echo 1 ;
    graph = 'N' ;

*** Paramètres ...

    rayon = 10. ;

    nb_div = 2 ;

    valyoun = 6.825e+7 ;
    valnu   = 0.3 ;
    lepaiss = 0.04 ;

*** Points ...

    dens (pi * rayon / 4 / nb_div) ;

    x1 = rayon ;
    x2 = (rayon * rayon / 2) ** 0.5 ;
    x3 = (rayon * rayon / 3) ** 0.5 ;

    p0 = 0. 0. 0. ;
    px = x1 0. 0. ;
    py = 0. x1 0. ;
    pz = 0. 0. x1 ;

    pxy = x2 x2 0. ;
    pxz = x2 0. x2 ;
    pyz = 0. x2 x2 ;

    pxyz = x3 x3 x3 ;

*** Lignes ...

    lixy = cerc nb_div px p0 pxy ;
    lixz = cerc nb_div px p0 pxz ;

    liyx = cerc nb_div py p0 pxy ;
    liyz = cerc nb_div py p0 pyz ;

    lizx = cerc nb_div pz p0 pxz ;
    lizy = cerc nb_div pz p0 pyz ;

    licxy = cerc nb_div pxyz p0 pxy ;
    licxz = cerc nb_div pxyz p0 pxz ;
    licyz = cerc nb_div pxyz p0 pyz ;

*** Surfaces ...

    opti elem qua4 ;

    quartx = dall lixy (inve licxy) licxz (inve lixz) sphe p0 ;
    quarty = dall liyz (inve licyz) licxy (inve liyx) sphe p0 ;
    quartz = dall lizx (inve licxz) licyz (inve lizy) sphe p0 ;

    su1 = quartx et quarty et quartz ;

    si (ega 'O' graph) ;
       titr 'Le maillage (un quart de la hemisphere)' ;
       trac su1 ;
    finsi ;

*** Modèle ...

    mo1 = mode su1 mecanique elastique coq4 ;
    ma1 = mate mo1 YOUN valyoun NU valnu EPAI lepaiss ;

*** Rigidité et CL ...

    ri1 = rigi mo1 ma1 ;
    cl1 = symt depl rota px p0 pz su1 ;
    cl2 = symt depl rota py p0 pz su1 ;
    cl3 = bloq uz pz ;

    rielas = ri1 et cl1 et cl2 et cl3 ;

    si(ega 'O'  graph) ;
       titr 'La structure avec ses conditions aux limites' ;
       trac rielas ;
    finsi ;

*** Chargement ...

    val_forc = 1. ;

    fo1 = forc fx     val_forc  px ;
    fo2 = forc fy (-1*val_forc) py ;

    fortot = fo1 et fo2 ;

*** Résolution ...

    de1 = reso rielas fortot ;

    si( ega 'O' graph) ;
       titr 'La structure deformee' ;
       vf1 = vect fortot FX FY FZ (rayon / val_forc / 5.) JAUN ;
       defo0 = defo su1 de1 0.     VERT ;
       defo1 = defo su1 de1    vf1 ROUG ;
       trac (defo0 et defo1) ;
    finsi ;

*** Test ...

    val_ref = 0.0925 ;

    dx =       extr de1 ux px  ;
    dy = -1 * (extr de1 uy py) ;

    errx = 100 * (dx - val_ref) / val_ref ;
    erry = 100 * (dy - val_ref) / val_ref ;

    opti echo 0 ;
    mess 'Déplacement de référence = ' val_ref ;
    mess 'Déplacement en X =' dx ' soit' errx ' % d erreur' ;
    mess 'Déplacement en Y =' dy ' soit' erry ' % d erreur' ;
    opti echo 1 ;

*** Bye ...

    fin ;