Télécharger la plaquette Cast3M


fr en

Télécharger frene0.eso

Retour à la liste

Numérotation des lignes : 

frene0

  1. C FRENE0    SOURCE    CB215821  23/11/02    21:15:04     11779          
  2.  
  3.       SUBROUTINE FRENE0(arc, repere)
  4. C ====================================================================
  5. C Sous-programme FRENE0 : calcul du repère de Frenet d'une ligne
  6. C
  7. C --------------------------
  8. C Paramètres Entrée/Sortie :
  9. C --------------------------
  10. C
  11. C E/  arc    : Ligne le long de laquelle on veut calculer le repère
  12. C
  13. C  /S repere : CHPOINT contenant les vecteurs du repère
  14. C ====================================================================
  15.  
  16.       IMPLICIT INTEGER(I-N)
  17.       IMPLICIT REAL*8(A-H,O-Z)
  18.  
  19. C --------------------------------------------------------------------
  20. C     Declarations
  21. C --------------------------------------------------------------------
  22.  
  23. -INC PPARAM
  24. -INC CCOPTIO
  25. -INC SMCOORD
  26. -INC SMELEME
  27. -INC SMCHPOI
  28. -INC CCREEL
  29.  
  30.       character*4 les_composantes(9)
  31.  
  32.       data les_composantes/'TX', 'TY', 'TZ', 'NX', 'NY', 'NZ',
  33.      $                     'BX', 'BY', 'BZ'/
  34.  
  35.       integer element, premier, dernier, composante_connexe, ecart
  36.  
  37.       real*8 oa(3), ob(3), oc(3),
  38.      $       ab(3), bc(3), axe(3), test(3),
  39.      $       tangent(3), normal(3), binormal(3),
  40.      $       norme_ab, norme_bc, norme_bin
  41.  
  42.       logical est_seg2, est_seg3, est_ouvert, est_transition,
  43.      $        sont_colineaires
  44.  
  45.       segment nombre_occurences(nbpts)
  46.       segment numero_du_poi1(nbpts)
  47.       segment le_premier_element(0)
  48.       segment le_dernier_element(0)
  49.  
  50.       pointeur repere.mchpoi
  51.       pointeur arc.meleme, arc_seg2.meleme, poi1_de_arc.meleme
  52.  
  53.  
  54. C --------------------------------------------------------------------
  55. C     Initialisations et preparation du travail
  56. C --------------------------------------------------------------------
  57.  
  58.       poi1_de_arc = arc
  59.       call change(poi1_de_arc, 1)
  60.       if(ierr.ne.0) return
  61.  
  62.       est_seg2 = (arc.itypel).eq.2
  63.       est_seg3 = (arc.itypel).eq.3
  64.       if(.not.(est_seg2.or.est_seg3)) then
  65.         call erreur(16)
  66.       endif
  67.  
  68.       arc_seg2 = arc
  69.       if(est_seg3) then
  70.         call change(arc_seg2, 2)
  71.         if(ierr.ne.0) return
  72.       endif
  73.  
  74. C      On verifie l'orientation des elements avec versen
  75. C      Pour ca on est oblige de faire des ecrobj/lirobj
  76. C      => Peut-etre faudrait-il creer un deuxieme niveau dans versen
  77. C      pour pouvoir l'appeler directement ?
  78.       call ecrobj('MAILLAGE', arc_seg2)
  79.       call versen
  80.       call lirobj('MAILLAGE', arc_seg2, 1, iretou)
  81.       if(ierr.ne.0) return
  82.  
  83. C     On verifie que le contour est simple
  84. C     (un noeud appartient a au plus deux elements)
  85. C     On fait comme dans ligmai.eso
  86.       segini nombre_occurences
  87.       maxel = 0
  88.       nombre_de_seg2 = arc_seg2.num(/2)
  89.       do element=1,nombre_de_seg2
  90.         do noeud_du_seg=1,2
  91.           noeud = arc_seg2.num(noeud_du_seg, element)
  92.           nombre_occurences(noeud) = nombre_occurences(noeud) + 1
  93.           maxel = max(maxel, nombre_occurences(noeud))
  94.         enddo
  95.       enddo
  96.       segsup nombre_occurences
  97.       if(maxel.gt.2) call erreur(426)
  98.  
  99. C     Initialisation du segment mpoval
  100.       nombre_noeuds = poi1_de_arc.num(/2)
  101.       n = nombre_noeuds
  102.       nombre_composantes = idim**2
  103.       nc = nombre_composantes
  104.       segini mpoval
  105.  
  106. C     Initialisation et remplissage de numero_du_poi1 (numerotation globale/locale)
  107.       segini numero_du_poi1
  108.       do i=1,nombre_noeuds
  109.         numero_du_poi1(poi1_de_arc.num(1, i)) = i
  110.       enddo
  111.       segdes poi1_de_arc
  112.  
  113. C     On compte le nombre de composantes connexes et on stocke
  114. C     les premier et dernier elements de chacune
  115.       segini le_premier_element, le_dernier_element
  116.       le_premier_element(**) = 1
  117.       do element=1,(nombre_de_seg2 - 1)
  118.         noeud_1 = arc_seg2.num(2, element)
  119.         noeud_2 = arc_seg2.num(1, element + 1)
  120.         est_transition = noeud_1.ne.noeud_2
  121.         if(est_transition) then
  122.           le_dernier_element(**) = element
  123.           le_premier_element(**) = element + 1
  124.         endif
  125.       enddo
  126.       le_dernier_element(**) = arc_seg2.num(/2)
  127.       nombre_composantes_connexes = le_premier_element(/1)
  128.  
  129.  
  130. C --------------------------------------------------------------------
  131. C     On fait le travail
  132. C --------------------------------------------------------------------
  133.  
  134. C     Comme un noeud a ses coordonnees + sa densite il faut prendre
  135. C     dim + 1 pour obtenir ses coordonnees a partir de xcoor
  136.       idim1 = idim + 1
  137.  
  138.       tangent = 0.D0
  139.       normal = 0.D0
  140.       binormal = (/ 0.D0, 0.D0, 1.D0 /)
  141.       ab = 0.D0
  142.       bc = 0.D0
  143.  
  144.       SEGACT,MCOORD
  145.       do composante_connexe=1,nombre_composantes_connexes
  146.  
  147. C     Travail sur tous les noeuds interieurs
  148. C     **************************************
  149.  
  150.         premier = le_premier_element(composante_connexe)
  151.         dernier = le_dernier_element(composante_connexe)
  152.  
  153.         nombre_iterations = dernier - premier + 1
  154.  
  155.         noeud_debut = arc_seg2.num(1, premier)
  156.         noeud_fin = arc_seg2.num(2, dernier)
  157.         est_ouvert = noeud_debut.ne.noeud_fin
  158.         if(est_ouvert) nombre_iterations = nombre_iterations - 1
  159.  
  160.         do element=premier,(premier + nombre_iterations - 1)
  161. C         Pour chaque 2eme noeud d'element on calcule d'abord le
  162. C         vecteur binormal, puis on estime le vecteur tangent,
  163. C         puis on calcule le vecteur normal par produit vectoriel,
  164. C         et enfin on recalcule le vecteur tangent par produit
  165. C         vectoriel pour etre sur que les 3 vecteurs sont orthogonaux
  166.  
  167.           ia = arc_seg2.num(1, element)
  168.           ib = arc_seg2.num(2, element)
  169.           if(.not.est_ouvert.and.element.eq.dernier) then
  170.             ic = arc_seg2.num(2, premier)
  171.           else
  172.             ic = arc_seg2.num(2, element + 1)
  173.           endif
  174.  
  175.           do i=1,idim
  176.             oa(i) = xcoor((ia-1)*idim1 + i)
  177.             ob(i) = xcoor((ib-1)*idim1 + i)
  178.             oc(i) = xcoor((ic-1)*idim1 + i)
  179.           enddo
  180.  
  181.           if(idim.eq.3) then
  182.             do i=1,3
  183.               ab(i) = ob(i) - oa(i)
  184.               bc(i) = oc(i) - ob(i)
  185.             enddo
  186.             call provec(ab, bc, binormal)
  187.  
  188.             call norme(ab, norme_ab)
  189.             call norme(bc, norme_bc)
  190.             call norme(binormal, norme_bin)
  191.             norme_bin = norme_bin + xpetit
  192.             sont_colineaires = norme_bin
  193.      $                              .le.(norme_ab*norme_bc*xzprec)
  194.             if(sont_colineaires) then
  195. C             On tente de prendre le vecteur binormal en A
  196.               do i=1,3
  197.                 binormal(i) = vpocha(numero_du_poi1(ia), 6 + i)
  198.               enddo
  199.               call norme(binormal, norme_bin)
  200.               norme_bin = norme_bin + xpetit
  201.               if(.not.a_egale_b(norme_bin, 1.D0)) then
  202. C               Il n'y a pas de vecteur unitaire en A donc on appelle perpen
  203.                 call perpen(ab, binormal)
  204.               endif
  205.             else
  206.               call normer(binormal)
  207.             endif
  208.           endif
  209.  
  210.           do i=1,idim
  211.             tangent(i) = oc(i) - oa(i)
  212.           enddo
  213.           call provec(binormal, tangent, normal)
  214.           call normer(normal)
  215.  
  216.           call provec(normal, binormal, tangent)
  217.  
  218.           noeud = numero_du_poi1(ib)
  219.           do i=1,idim
  220.             vpocha(noeud, i) = tangent(i)
  221.             vpocha(noeud, idim + i) = normal(i)
  222.             if(idim.eq.3) vpocha(noeud, 6 + i) = binormal(i)
  223.           enddo
  224.  
  225.         enddo
  226.  
  227. C     Travail sur les extermites (si le contour est ouvert)
  228. C     ****************************************************
  229.  
  230.         if(est_ouvert) then
  231. C         On fait une boucle pour traiter les deux extremites :
  232. C         - on calcule le vecteur tangent
  233. C         - on calcule la binormale si on est en 3D (par orthonormalisation
  234. C           de la valeur voisine ou en appelant perpen si contour = 1 element)
  235. C         - on calcule la normale par le produit vectoriel des deux vecteurs precedents
  236. C           (attention a la boucle, si contour = 1 element alors dernier - premier = 0 !)
  237.           ecart = max(1, dernier - premier)
  238.           do element=premier,dernier,ecart
  239.             ia = arc_seg2.num(1, element)
  240.             ib = arc_seg2.num(2, element)
  241.             if(element.eq.premier) then
  242.               noeud_calcul = ia
  243.               noeud_voisin = ib
  244.             else
  245.               noeud_calcul = ib
  246.               noeud_voisin = ia
  247.             endif
  248.  
  249.             do i=1,idim
  250.               oa(i) = xcoor((ia-1)*idim1 + i)
  251.               ob(i) = xcoor((ib-1)*idim1 + i)
  252.               tangent(i) = ob(i) - oa(i)
  253.             enddo
  254.             call normer(tangent)
  255.  
  256.             if(idim.eq.3) then
  257.               if(premier.ne.dernier) then
  258.                 noeud = numero_du_poi1(noeud_voisin)
  259.                 do i=1,3
  260.                   binormal(i) = vpocha(noeud, 6 + i)
  261.                 enddo
  262.                 call scal(tangent, binormal, t_scal_b)
  263.                 do i=1,3
  264.                   binormal(i) = binormal(i) - t_scal_b * tangent(i)
  265.                 enddo
  266.                 call normer(binormal)
  267.               else
  268.                 call perpen(tangent, binormal)
  269.               endif
  270.             endif
  271.  
  272.             call provec(binormal, tangent, normal)
  273.  
  274.             noeud = numero_du_poi1(noeud_calcul)
  275.             do i=1,idim
  276.               vpocha(noeud, i) = tangent(i)
  277.               vpocha(noeud, idim + i) = normal(i)
  278.               if(idim.eq.3) vpocha(noeud, 6 + i) = binormal(i)
  279.             enddo
  280.           enddo
  281.  
  282. C         Si un seul element, on met les memes valeurs au deuxieme noeud
  283.           if(premier.eq.dernier) then
  284.             noeud = numero_du_poi1(noeud_voisin)
  285.             do i=1,idim
  286.               vpocha(noeud, i) = tangent(i)
  287.               vpocha(noeud, idim + i) = normal(i)
  288.               if(idim.eq.3) vpocha(noeud, 6 + i) = binormal(i)
  289.             enddo
  290.           endif
  291.         endif
  292.  
  293.       enddo
  294.       SEGDES,MCOORD
  295.  
  296.       segsup numero_du_poi1, arc_seg2
  297.       segsup le_premier_element, le_dernier_element
  298.  
  299. C     On finalise les segments du CHPOINT
  300.       nsoupo = 1
  301.       nat = 1
  302.       segini repere
  303.       repere.jattri(1) = 1
  304.       repere.ifopoi = ifour
  305.       repere.mtypoi = ' '
  306.       repere.mochde = 'REPERE DE FRENET'
  307.       segini msoupo
  308.       repere.ipchp(1) = msoupo
  309.       igeoc = poi1_de_arc
  310.       do j=1,idim
  311.         do i=1,idim
  312.           nocomp(idim*(j-1) + i) = les_composantes(3*(j-1) + i)
  313.           noharm(idim*(j-1) + i) = nifour
  314.         enddo
  315.       enddo
  316.       ipoval = mpoval
  317.       segdes repere, msoupo, mpoval
  318.  
  319.       END
  320.  
  321.  
  322.  
  323.  
  324.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales