Esope source files

   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 :