Télécharger la plaquette Cast3M


fr en

Télécharger rlexvb.eso

Retour à la liste

Numérotation des lignes : 

rlexvb

  1. C RLEXVB    SOURCE    OF166741  24/12/13    21:17:33     12097          
  2.       SUBROUTINE RLEXVB(MELFL,MELFP,MELSOM,MLELEM)
  3. C
  4.       IMPLICIT INTEGER(I-N)
  5.       INTEGER  NSOMM, LAST, NTOTN, NTOTE, NBSOUS, NBELEM, NBNO
  6.      &       , IELEM, NGF, NGF1, INOEU, NGS1, NLS1, ISOUS
  7.      &       , IPOS1, IPOSV1, NVOIS1, INOEU2, NGS2, NLS2
  8.      &       , IVOIS, IPOS, IELEMF, I1
  9.       LOGICAL LOGCON
  10. C
  11. -INC SMCOORD
  12. -INC SMELEME
  13.       INTEGER JG
  14. -INC SMLENTI
  15.  
  16. -INC PPARAM
  17. -INC CCOPTIO
  18. C
  19.       INTEGER NBL, NBTPOI
  20.       SEGMENT MLELEM
  21.       INTEGER INDEX(NBL+1)
  22.       INTEGER LESPOI(NBTPOI)
  23.       ENDSEGMENT
  24. C
  25.       POINTEUR MELSOM.MELEME, MELFL.MELEME, MELFP.MELEME
  26.      &         ,MLESOM.MLENTI, MLEFP.MLENTI, MPOS.MLENTI, MTOUC.MLENTI
  27.      &         ,MLELE1.MLELEM, MELFP1.MELEME
  28. C
  29. C**** Le MELEME SOMMET
  30. C
  31.       CALL KRIPAD(MELSOM,MLESOM)
  32. C
  33. C     MLESOM: numerotation globale -> locale
  34. C
  35. C**** En KRIPAD
  36. C     SEGACT MELSOM
  37. C     SEGINI MLESOM
  38. C
  39.       NSOMM = MELSOM.NUM(/2)
  40.       JG=NSOMM
  41.       SEGINI MTOUC
  42. C     MTOUC.LECT(NLS1) = surestimation de nombre de voisins de NLS1
  43.       SEGINI MPOS
  44. C     MPOS.LECT + LAST : liste chaînée des sommets sur le bord
  45.       LAST=-1
  46. C
  47.       SEGACT MELFP
  48. C
  49. C**** En 2D MELFP est un maillage elementaire
  50. C     En 3D pas à priori
  51. C
  52. C     NTOTE : nombre total d'elts dans la même structure
  53. C
  54.       NTOTN = 0
  55.       NTOTE = 0
  56.       NBSOUS=MELFP.LISOUS(/1)
  57. C     NBSOUS=0 fais un peux chier!
  58.       JG=MAX(NBSOUS,1)
  59.       SEGINI MLEFP
  60.       IF(NBSOUS .EQ. 0)THEN
  61.          MLEFP.LECT(1)=MELFP
  62.       ELSE
  63.          DO ISOUS=1,NBSOUS,1
  64.             MLEFP.LECT(ISOUS)=MELFP.LISOUS(ISOUS)
  65.          ENDDO
  66.       ENDIF
  67.       SEGDES MELFP
  68. C
  69.       SEGACT MELFL
  70.       NBSOUS=MELFL.LISOUS(/1)
  71.       IF(NBSOUS .NE. 0)THEN
  72.          WRITE(IOIMP,*) 'FACEL = ???'
  73.          WRITE(IOIMP,*) 'Subroutine rlexvb.eso'
  74.          CALL ERREUR(5)
  75.          GOTO 9999
  76.       ENDIF
  77. C
  78.       NBSOUS=JG
  79.       IELEMF=0
  80.       DO ISOUS = 1, NBSOUS, 1
  81.          MELFP1=MLEFP.LECT(ISOUS)
  82.          SEGACT MELFP1
  83.          NBELEM=MELFP1.NUM(/2)
  84.          NBNO=MELFP1.NUM(/1) - 1
  85.          DO IELEM = 1, NBELEM,1
  86.             IELEMF=IELEMF+1
  87.             NGF=MELFP1.NUM(NBNO+1,IELEM)
  88.             NGF1=MELFL.NUM(2,IELEMF)
  89.             IF(NGF .NE. NGF1)THEN
  90.                WRITE(IOIMP,*) 'FACEL = ???'
  91.                WRITE(IOIMP,*) 'Subroutine rlexvb.eso'
  92.                CALL ERREUR(5)
  93.                GOTO 9999
  94.             ENDIF
  95.             IF(MELFL.NUM(1,IELEMF) .EQ. MELFL.NUM(3,IELEMF))THEN
  96. C
  97. C************* La face est sur le bord
  98. C
  99.                DO INOEU = 1, NBNO, 1
  100.                   NGS1 = MELFP1.NUM(INOEU,IELEM)
  101.                   NLS1 = MLESOM.LECT(NGS1)
  102.                   IF(MPOS.LECT(NLS1).EQ.0)THEN
  103.                      NTOTE=NTOTE+1
  104.                      MPOS.LECT(NLS1)=LAST
  105.                      LAST=NLS1
  106.                   ENDIF
  107.                ENDDO
  108.             ENDIF
  109.          ENDDO
  110.       ENDDO
  111. C
  112. C     NTOTN : nombre total de noeuds dans la structure sommets de bords
  113. C             + sommets voisins
  114. C     N.B. Les sommets voisins peuvent ne pas etre sur le bord
  115. C
  116.       NTOTN = 0
  117.       DO ISOUS = 1, NBSOUS, 1
  118.          MELFP1=MLEFP.LECT(ISOUS)
  119.          NBELEM=MELFP1.NUM(/2)
  120.          NBNO=MELFP1.NUM(/1) - 1
  121.          DO IELEM = 1, NBELEM,1
  122.             DO INOEU = 1, NBNO, 1
  123.                NGS1 = MELFP1.NUM(INOEU,IELEM)
  124.                NLS1 = MLESOM.LECT(NGS1)
  125.                IF(MPOS.LECT(NLS1).NE.0)THEN
  126. C
  127. C******************** Le sommet est sur le bord
  128. C
  129.                   MTOUC.LECT(NLS1)=MTOUC.LECT(NLS1)+NBNO -1
  130.                   NTOTN=NTOTN+(NBNO - 1)
  131.                ENDIF
  132.             ENDDO
  133.          ENDDO
  134.       ENDDO
  135. C
  136. C**** On stoke l'informaton dans une LISTE SEQUENTIELLE INDEXEE D'ELEMENTS
  137. C
  138. C     NBL      : NOMBRE D'ELEMENTS
  139. C     NBTPOI   : NOMBRE TOTAL DE POINTS REFERENCEES
  140. C     INDEX(I) : INDICE DU 1ER POINT DU IEME ELEMENT
  141. C                DANS LE TABLEAU LESPOI
  142. C     LESPOI(INDEX(I) -> INDEX(I+1)-1) : NUMERO DES NOEUDS
  143. C                DU IEME ELEMENT
  144. C
  145. C     En plus, après le REPEAT UNTIL,
  146. C     MPOS.LECT(NLS1) = position du sommet NLS1 i.e.
  147. C        J=MPOS.LECT(NLS1) -> NLS1=LESPOI(INDEX(J))
  148. C                            se voisins sont dedans LESPOI(INDEX(J)+1 -> INDEX(I+1)-1)
  149.       NBL=NTOTE
  150.       NBTPOI=NBL+NTOTN
  151.       SEGINI MLELEM
  152.       IPOS1=0
  153.       MLELEM.INDEX(1)=1
  154. C
  155. C**** REPEAT UNTIL
  156. C
  157.  1    CONTINUE
  158.       NLS1 = LAST
  159.       IF(NLS1 .NE. -1)THEN
  160.          IPOS1=IPOS1+1
  161.          MLELEM.LESPOI(MLELEM.INDEX(IPOS1))=NLS1
  162.          LAST = MPOS.LECT(NLS1)
  163.          MPOS.LECT(NLS1)=IPOS1
  164.          MLELEM.INDEX(IPOS1+1)=MLELEM.INDEX(IPOS1)+MTOUC.LECT(NLS1)+1
  165.          MTOUC.LECT(NLS1)=0
  166.          GOTO 1
  167.       ENDIF
  168. C
  169. C**** N.B. MTOUC.LECT(I) = 0 \forall I
  170. C
  171. C     Fin REPEAT UNTIL
  172. C     On remplie LESPOI
  173. C
  174.       DO ISOUS = 1, NBSOUS, 1
  175.          MELFP1 = MLEFP.LECT(ISOUS)
  176.          NBELEM=MELFP1.NUM(/2)
  177.          NBNO=MELFP1.NUM(/1) - 1
  178.          DO IELEM = 1, NBELEM, 1
  179.             DO INOEU = 1, NBNO, 1
  180.                NGS1 = MELFP1.NUM(INOEU,IELEM)
  181.                NLS1 = MLESOM.LECT(NGS1)
  182. C              IPOS1 = 0 -> NGS1 n'est pas sur le bords
  183. C              IPOS1 = C > 0
  184. C                C= position de NGS1 (NLS1) dedans MLELEM.INDEX
  185.                IPOS1 = MPOS.LECT(NLS1)
  186.                IF(IPOS1 .GT. 0)THEN
  187. C
  188. C                 IPOSV1 = position du premier voisin de NLS1 dedans
  189. C                       MLEMEM.LESPOI
  190.                   IPOSV1 = MLELEM.INDEX(IPOS1)+1
  191. C                 NVOIS1 = surestimation du nombre de voisins de NLS1
  192.                   NVOIS1 = MLELEM.INDEX(IPOS1+1)-IPOSV1
  193.                   DO INOEU2 = 1, NBNO, 1
  194.                      NGS2=MELFP1.NUM(INOEU2,IELEM)
  195.                      IF(NGS1 .NE. NGS2)THEN
  196.                         NLS2 = MLESOM.LECT(NGS2)
  197.                         IVOIS=0
  198. C
  199. C**************** REPEAT UNTIL
  200. C
  201.  10                     CONTINUE
  202.                         IF(MLELEM.LESPOI(IPOSV1+IVOIS) .EQ. 0)THEN
  203.                            MTOUC.LECT(NLS1)=MTOUC.LECT(NLS1)+1
  204.                            MLELEM.LESPOI(IPOSV1+IVOIS) = NLS2
  205.                         ELSE
  206.                            LOGCON = (NLS2 .NE. MLELEM.LESPOI(IPOSV1
  207.      $                          +IVOIS))
  208.                            IVOIS=IVOIS+1
  209.                            IF(LOGCON)THEN
  210.                               IF(IVOIS .GT. NVOIS1)THEN
  211. C
  212. C************************* Il y a un erreur de programmation (NVOIS1 n'a
  213. C                          pas bien été évalué)
  214. C
  215.                                  WRITE(IOIMP,*) 'Subroutine rlexvb.eso'
  216.                                  CALL ERREUR(5)
  217.                                  GOTO 9999
  218.                               ELSE
  219.                                  GOTO 10
  220.                               ENDIF
  221.                            ENDIF
  222.                         ENDIF
  223.                      ENDIF
  224.                   ENDDO
  225.                ENDIF
  226.             ENDDO
  227.          ENDDO
  228.          SEGDES MELFP1
  229.       ENDDO
  230.       SEGSUP MLEFP
  231. C
  232. C**** LESPOI SURDIMENSIONÉ
  233. C
  234.       NBTPOI=0
  235.       DO IELEM = 1 , NTOTE, 1
  236.          NLS1 = MLELEM.LESPOI(MLELEM.INDEX(IELEM))
  237.          NBTPOI=NBTPOI+1+MTOUC.LECT(NLS1)
  238.       ENDDO
  239.       NBL=NTOTE
  240.       SEGINI MLELE1
  241. C
  242.       MLELE1.INDEX(1)=1
  243.       DO IELEM = 1, NBL , 1
  244.          IPOS=MLELEM.INDEX(IELEM)
  245.          IPOS1=MLELE1.INDEX(IELEM)
  246.          NLS1 =  MLELEM.LESPOI(IPOS)
  247.          NGS1 =  MELSOM.NUM(1,NLS1)
  248.          NVOIS1 = MTOUC.LECT(NLS1)
  249.          MLELE1.LESPOI(IPOS1)=NGS1
  250.          MLELE1.INDEX(IELEM+1)=MLELE1.INDEX(IELEM)+NVOIS1+1
  251.          DO INOEU = 1, NVOIS1, 1
  252.             NLS2 = MLELEM.LESPOI(IPOS+INOEU)
  253.             IF(NLS2 .EQ.0)THEN
  254.                WRITE(IOIMP,*) 'Subroutine rlexvb.eso'
  255.                CALL ERREUR(5)
  256.                GOTO 9999
  257.             ELSE
  258.                NGS2=MELSOM.NUM(1,NLS2)
  259.                MLELE1.LESPOI(IPOS1+INOEU)=NGS2
  260.             ENDIF
  261.          ENDDO
  262.       ENDDO
  263. C
  264.       SEGSUP MLELEM
  265.       MLELEM=MLELE1
  266. C
  267.       SEGSUP MLESOM
  268.       SEGSUP MPOS
  269.       SEGSUP MTOUC
  270.       SEGDES MLELEM
  271.       SEGDES MELSOM
  272.       SEGDES MELFL
  273. C
  274. C**** Problème détecté dans un cas test: mal-conditionnent
  275. C     de la matrice associé à la méthode de moindres carres
  276. C     aux bords.
  277. C     Quand possible, on elimine de MLELEM tous les voisins
  278. C     sommets qui sont sur le bord -> MLELE1
  279. C
  280.       SEGACT MLELEM
  281.       NBL=MLELEM.INDEX(/1)-1
  282.       NBTPOI=MLELEM.LESPOI(/1)
  283.       SEGINI MLELE1
  284. C
  285. C**** MLESOM.MLENTI contient les sommets aux bords
  286. C
  287.       JG=nbpts
  288.       SEGINI MLESOM
  289.       DO IELEM=1,NBL,1
  290.          IPOS=MLELEM.INDEX(IELEM)
  291.          NGS1=MLELEM.LESPOI(IPOS)
  292.          MLESOM.LECT(NGS1)=IELEM
  293.       ENDDO
  294. C
  295. C
  296. C**** On rempli MLELE1
  297. C
  298.       NBTPOI=0
  299.       NVOIS1=0
  300. C
  301. C**** Pour chaque SOMMET NGS1, NVOIS1 = nombre des voisins qui ne
  302. C     sont pas sur le bords
  303. C
  304.       IPOS1=MLELEM.INDEX(1)
  305.       DO IELEM=1,NBL,1
  306.          IPOS=IPOS1
  307.          IPOS1=MLELEM.INDEX(1+IELEM)
  308.          NBTPOI=NBTPOI+1
  309.          MLELE1.INDEX(IELEM)=NBTPOI
  310.          NGS1=MLELEM.LESPOI(IPOS)
  311.          MLELE1.LESPOI(NBTPOI)=NGS1
  312.          DO I1=IPOS+1,IPOS1-1,1
  313.             NGS2=MLELEM.LESPOI(I1)
  314.             NLS2=MLESOM.LECT(NGS2)
  315.             IF(NLS2.EQ.0)THEN
  316. C
  317. C************* NGS2 n'est pas sur le bord
  318. C
  319.                NVOIS1=NVOIS1+1
  320.                MLELE1.LESPOI(NBTPOI+NVOIS1)=NGS2
  321.             ENDIF
  322.          ENDDO
  323.          IF(NVOIS1.EQ.0)THEN
  324. C
  325. C********** Tous les voisins de NGS1 sont sur le bords
  326. C
  327.             DO I1=IPOS+1,IPOS1-1,1
  328.                NGS2=MLELEM.LESPOI(I1)
  329.                NVOIS1=NVOIS1+1
  330.                MLELE1.LESPOI(NBTPOI+NVOIS1)=NGS2
  331.             ENDDO
  332.          ENDIF
  333.          NBTPOI=NBTPOI+NVOIS1
  334.          NVOIS1=0
  335.       ENDDO
  336.       MLELE1.INDEX(NBL+1)=NBTPOI+1
  337. C
  338.       SEGSUP MLELEM
  339.       SEGADJ MLELE1
  340.       MLELEM=MLELE1
  341.       SEGDES MLELEM
  342.       SEGSUP MLESOM
  343. C
  344.  9999 RETURN
  345.       END
  346.  
  347.  
  348.  
  349.  
  350.  

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