Télécharger la plaquette Cast3M


fr en

Télécharger dualis.eso

Retour à la liste

Numérotation des lignes : 

dualis

  1. C DUALIS    SOURCE    CHAT      05/01/12    22:58:42     5004
  2. C   dualise un maillage pour le maillage par polygone
  3. C
  4.       subroutine dualis(fer,fer1,xpro,xproj1)
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8 (A-H,O-Z)
  7.       SEGMENT /FER/(NFI(ITT),MAI(IPP),ITOUR)
  8.       SEGMENT XPRO
  9.         REAL*8 XPROJ(3,1)
  10.       ENDSEGMENT
  11.       POINTEUR XPROJ1.XPRO, FER1.FER
  12.       xproj1=xpro
  13.       segini,xpro=xproj1
  14.       segini,fer1=fer
  15.       do 10 it=1,itour
  16.         ideb=mai(it)+1
  17.         ifin=mai(it+1)
  18.         do 20 ip1=ideb,ifin
  19.           ip2=ip1+1
  20.           if (ip2.gt.ifin) ip2=ideb
  21.           xproj(1,nfi(ip1))=(xproj1.XPROJ(1,nfi(ip1))+
  22.      >                       xproj1.XPROJ(1,nfi(ip2)))/2
  23.           xproj(2,nfi(ip1))=(xproj1.XPROJ(2,nfi(ip1))+
  24.      >                       xproj1.XPROJ(2,nfi(ip2)))/2
  25.           xproj(3,nfi(ip1))=(xproj1.XPROJ(3,nfi(ip1))+
  26.      >                       xproj1.XPROJ(3,nfi(ip2)))/2
  27.  
  28.   20    continue
  29.   10  continue
  30. *
  31. *     Decalage des points duals associes aux sommets
  32. *     interieurs d'un contour concave
  33. *
  34.       do 40 it=1,itour
  35.         ideb=mai(it)+1
  36.         ifin=mai(it+1)
  37.         do 30 ip1=ideb,ifin
  38.           ip2=ip1+1
  39.           if (ip2.gt.ifin) ip2=ideb
  40.           ip3 = ip2 + 1
  41.           if (ip3.gt.ifin) ip3=ideb
  42. *
  43. *         Produit vectoriel des cotes liant un sommet
  44. *
  45.           xvec1 = xproj1.XPROJ(1,nfi(ip2)) - xproj1.XPROJ(1,nfi(ip1))
  46.           yvec1 = xproj1.XPROJ(2,nfi(ip2)) - xproj1.XPROJ(2,nfi(ip1))
  47.           xvec2 = xproj1.XPROJ(1,nfi(ip3)) - xproj1.XPROJ(1,nfi(ip2))
  48.           yvec2 = xproj1.XPROJ(2,nfi(ip3)) - xproj1.XPROJ(2,nfi(ip2))
  49.  
  50.           prod = xvec1 * yvec2 - xvec2 * yvec1
  51.           IF ( prod .LT. -1.E-10) THEN
  52. *
  53. *           Le sommet est "entrant", decalage des points duals voisins.
  54. *
  55.             xproj(1,nfi(ip1))= xproj(1,nfi(ip1)) + (xvec1 - xvec2)/4
  56.             xproj(2,nfi(ip1))= xproj(2,nfi(ip1)) + (yvec1 - yvec2)/4
  57.             xproj(1,nfi(ip2))= xproj(1,nfi(ip2)) + (xvec1 - xvec2)/4
  58.             xproj(2,nfi(ip2))= xproj(2,nfi(ip2)) + (yvec1 - yvec2)/4
  59.  
  60.           ENDIF
  61.   30    continue
  62.   40  continue
  63.       end
  64.  
  65.  
  66.  

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