Télécharger la plaquette Cast3M


fr en

Télécharger mfilte.eso

Retour à la liste

Numérotation des lignes : 

mfilte

  1. C MFILTE    SOURCE    FD218221  25/12/15    21:15:02     12413          
  2.  
  3. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  4. CC Subroutine for generating a (rigidity) matrix used for filtering fields
  5. CC
  6. CC Authors:
  7. CC Guenhael Le Quilliec (LaMe - Polytech Tours)
  8. CC Thomas Fournier (Internship 2020 at LaMe - Polytech Tours under the supervision of G. Le Quilliec)
  9. CC
  10. CC Version:
  11. CC 1.0 2021/01/04 Original version (based on previous source file matrigv5.eso)
  12. CC 1.1 2021/01/07 The size of the cells is not anymore an option
  13. CC 2.0 2025/12/03 The matrix is no longer normalized (by dividing by the sum of the row)
  14. CC                This operation can now be performed using the * and NFIL operator
  15. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  16.  
  17. C     TODO GLQ : pour répartir plus rapidement les noeuds dans les cellules,
  18. C          retirer Xmin, Ymin et Zmin aux coordonnées X, Y et Z de
  19. C          tous les noeuds puis diviser par R et arrondir.
  20. C          On obtient ainsi les indices des cellules de chaque noeud.
  21. C          Sans doute plus rapide que l'approche actuelle?
  22.  
  23. C     TODO GLQ : tirer profit de la parallélisation pour accélérer la
  24. C          construction de la matrice.
  25.  
  26. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  27. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  28. CC
  29. CC  Initialisation de la subroutine
  30. CC
  31. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  32. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  33.  
  34.       SUBROUTINE MFILTE
  35.  
  36. C     Ces lignes indique que toutes les variables dont le nom commence par
  37. C     une lettre comprise entre I et N soient des entiers et
  38. C     celles avec un nom commençant par une lettre entre A et H ou O et Z
  39. C     soient des flottants (sur 8 octets / double précision)
  40.       IMPLICIT INTEGER  (I-N)
  41.       IMPLICIT REAL*8(A-H,O-Z)
  42.  
  43. C     Inclusion des segments standards utilisés dans cette subroutine
  44. C     ***************************************************************
  45.  
  46. -INC PPARAM
  47. C  Segment des options de calcul (IFOUR, etc.)
  48. -INC CCOPTIO
  49. C  Segment maillage
  50. -INC SMELEME
  51. C  Segment champ par point
  52. -INC SMCHPOI
  53. C  Segment liste d'entiers
  54. -INC SMLENTI
  55. C  Segment liste de réels
  56. -INC SMLREEL
  57. C  Segment coordonnées des noeuds
  58. -INC SMCOORD
  59. C  Segment matrice de rigidité
  60. -INC SMRIGID
  61.  
  62. C     Déclaration de segments complémentaires dédiés à cette subroutine
  63. C     *****************************************************************
  64.  
  65. C     Grille (type tableau 2D ou 3D) contenant les pointeurs des listes
  66. C     de noeuds du maillage contenus dans chacune des cellules.
  67. C     Comme la taille de grille n'est pas connue à l'avance,
  68. C     on passe par la définition d'un segment pour limiter
  69. C     l'occupation mémoire.
  70.       SEGMENT GRID2SEG
  71.           INTEGER GRID2(NX0,NY0)
  72.       ENDSEGMENT
  73.       SEGMENT GRID3SEG
  74.           INTEGER GRID3(NX0,NY0,NZ0)
  75.       ENDSEGMENT
  76.  
  77. C     Tableau contenant les indices maximum suivant X des cellules
  78. C     voisines dans un rayon R0 d'une cellule de référence
  79. C     Comme la taille de ce tableau n'est pas connue à l'avance,
  80. C     on passe par la définition d'un segment pour limiter
  81. C     l'occupation mémoire.
  82.       POINTEUR INDEX2.MLENTI
  83.       SEGMENT INDEX3SEG
  84.           INTEGER INDEX3(JG,JG)
  85.       ENDSEGMENT
  86.  
  87. C     Diverses variables
  88. C     ******************
  89.  
  90. C     Noms des inconnues primales et duales
  91.       CHARACTER*8 NOMP0, NOMD0
  92.  
  93. C     Diverses initialisations
  94. C     ************************
  95.  
  96. C     On indique une fois pour toutes le nombre d'inconnues duales
  97. C     des sous-matrices de rigidité qui seront créées (toujours 1)
  98.       NLIGRD = 1
  99.  
  100. C     Initialisation du nombre  de sous-matrices
  101.       NRIGEL = 0
  102.  
  103. C     Initialisation du nombre maximum de voisins
  104.       NBVMAX0 = 0
  105.  
  106. C     Initialisation du nombre de cellules non vides
  107.       NBNEC0 = 0
  108.  
  109. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  110. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  111. CC
  112. CC  Lecture des données d'entrée, traitement et vérifications
  113. CC
  114. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  115. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  116.  
  117. C     Lecture des données d'entrée
  118. C     ****************************
  119.  
  120. C     Lecture du champ par point des volumes
  121.       CALL LIROBJ('CHPOINT ',MCHPOI,1,IRETOU)
  122.       CALL ACTOBJ('CHPOINT ',MCHPOI,1)
  123.       IF (IERR.NE.0) RETURN
  124.  
  125. C     Lecture des données optionnelles :
  126. C     - le rayon R
  127.       R0 = 0.0
  128.       CALL LIRREE(R0,0,IRETOU)
  129.       IF (IERR.NE.0) RETURN
  130. C     - l'exposant MU
  131.       EMU0 = 1.0
  132.       CALL LIRREE(EMU0,0,IRETOU)
  133.       IF (IERR.NE.0) RETURN
  134. C     - la borne minimale des poids à ignorer TODO effectuer une vérification que sa valeur est bien entre 0.0 et 1.0
  135.       CRI0 = 0.0
  136.       CALL LIRREE(CRI0,0,IRETOU)
  137.       IF (IERR.NE.0) RETURN
  138. C     - la taille d'une cellule pour le classement des noeuds
  139.       C0 = R0
  140. c      CALL LIRREE(C0,0,IRETOU) >>> On supprime cette option car C0=R0 donne les meilleurs résultats dans la plupart des cas
  141. C     - le nom des inconnues primales
  142.       NOMP0 = 'SCAL'
  143.       CALL LIRCHA(NOMP0,0,IRETOU)
  144.       IF (IERR.NE.0) RETURN
  145. C     - le nom des inconnues duales
  146.       NOMD0 = NOMP0
  147.       CALL LIRCHA(NOMD0,0,IRETOU)
  148.       IF (IERR.NE.0) RETURN
  149.  
  150. C     traitement et vérifications
  151. C     ***************************
  152.  
  153. C     Activation du segment MCOORD de coordonnées des noeuds XCOOR
  154.       SEGACT MCOORD
  155.  
  156. C     On vérifie qu'il contient bien une et une seule partition géométrique
  157.       IF(IPCHP(/1).NE.1) THEN
  158. C         Si le champ est vide, on renvoie un entier égal à 1
  159.           IF(IPCHP(/1).EQ.0) THEN
  160.               CALL ECRREE(1.0d0)
  161.               RETURN
  162.           ENDIF
  163. C         Si il y a plus d'une partition on génère une erreur
  164. C         "Donnees incompatibles"
  165.           CALL ERREUR(21)
  166.       ENDIF
  167.  
  168. C     On suppose qu'il ne repose que sur une partition géométrique unique
  169. C     dont on récupère le sous-champ TODO à améliorer pour le cas où il y en a plusieurs
  170.       MSOUPO = IPCHP(1)
  171. C     On vérifie qu'il n'y a bien qu'une composante sinon on génère une erreur
  172. C     "Il faut specifier un champ par point avec une seule composante"
  173.       IF (NOCOMP(/2).NE.1) CALL ERREUR(180)
  174. C     On récupère le pointeur du tableau de valeurs du sous-champ
  175.       MPOVAL = IPOVAL
  176. C     On récupère le pointeur du maillage de la partition géométrique
  177.       MELEME = IGEOC
  178. C     On récupère le nombre de noeuds du maillage
  179.       NBN0 = NUM(/2)
  180.  
  181. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  182. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  183. CC
  184. CC  Si le rayon de filtrage est nul,
  185. CC  création d'une matrice de rigidité 'identité'
  186. CC
  187. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  188. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  189.  
  190.       IF(R0.GT.0.0) GO TO 2
  191.  
  192. C     La matrice contiendra une seule sous-matrice de NBN0 matrices élémentaires
  193.       NRIGEL = 1
  194.       NELRIG = NBN0
  195. C     On indique son nombre d'inconnues duales (déjà fait au début)
  196. C      NLIGRD = 1
  197. C     Elle possède une inconnue primale
  198.       NLIGRP = 1
  199. C     Le maillage correspondant aura NBN0 éléments à 1 noeud
  200.       NBNN = 1
  201.       NBELEM = NBN0
  202.       NBSOUS = 0
  203.       NBREF = 0
  204. C     Initialisation des segments matrice de rigidité, sous-matrice, descripteur et maillage
  205.       SEGINI MRIGID, XMATRI, DESCR, IPT1
  206. C     On indique son type de matrice
  207.       MTYMAT = 'RIGIDITE'
  208. C     On applique le mode de calcul du CHPOINT fourni
  209.       IFORIG = mchpoi.IFOPOI
  210. C     L'inconnue duale et l'inconnue primale porteront sur le premier (et seul)
  211. C     noeud de l'élément correspondant
  212.       NOELED(1) = 1
  213.       NOELEP(1) = 1
  214. C     Nom des inconnues
  215.       LISINC(1) = NOMP0
  216.       LISDUA(1) = NOMD0
  217. C     On indique que la matrice ne possède pas de symétries
  218.       SYMRE = 2
  219. C     On met à 1 les poids contenus dans la sous-matrice
  220.       DO 1 I = 1, NBN0
  221.           IPT1.NUM(1,I) = NUM(1,I)
  222.           RE(1,1,I) = 1.0
  223. 1     CONTINUE
  224. C     Affectation du maillage, du descripteur et du contenu de la sous-matrice
  225.       IRIGEL(1,1) = IPT1
  226.       IRIGEL(3,1) = DESCR
  227.       IRIGEL(4,1) = XMATRI
  228. C     On indique que la matrice globale ne possède pas de symétries
  229.       IRIGEL(7,1) = 2
  230. C     On précise son coefficient multiplicateur
  231.       COERIG(1) = 1
  232.  
  233. C     Passage direct à la fin de la subroutine
  234.       GO TO 99
  235.  
  236. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  237. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  238. CC
  239. CC  Sinon, création d'une matrice de rigidité
  240. CC  dont les sous-matrices contiennent chacune toutes
  241. CC  les matrices élémentaires de mêmes tailles
  242. CC  (nombre de noeuds influents identiques)
  243. CC
  244. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  245. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  246.  
  247. 2     CONTINUE
  248.  
  249. C     Initialisation de diverses listes
  250. C     *********************************
  251.  
  252. C     Liste d'entiers contenant les numéros des noeuds d'une cellule
  253.       POINTEUR MLENTC.MLENTI
  254.  
  255. C     Liste d'entiers contenant les numéros des noeuds du voisinage d'une cellule
  256.       POINTEUR MLENTV.MLENTI
  257.  
  258. C     Liste d'entiers contenant les numéros des voisins du noeud en cours de traitement
  259.       POINTEUR MLENTN.MLENTI
  260.  
  261. C     Liste de réels contenant les poids attribués aux voisins du noeud en cours de traitement
  262.       POINTEUR MLREEP.MLREEL
  263.  
  264. C     Liste des indices des sous-matrices de rigidité non vides
  265.       POINTEUR MLENTS.MLENTI
  266. C     Initialisation du segment
  267.       JG = 0
  268.       SEGINI MLENTS
  269.  
  270. C     Listes d'entiers pour stocker les matrices élémentaires, maillages et descripteurs des différentes sous-matrices
  271.       POINTEUR MLENTX.MLENTI
  272.       POINTEUR MLENTM.MLENTI
  273.       POINTEUR MLENTD.MLENTI
  274. C     Initialisation des segments
  275.       JG = 0
  276.       SEGINI MLENTX, MLENTM, MLENTD
  277.  
  278. C     Listes des indices des cellules non vides
  279.       POINTEUR NECX.MLENTI
  280.       POINTEUR NECY.MLENTI
  281.       POINTEUR NECZ.MLENTI
  282.       JG = 0
  283.       SEGINI NECX, NECY, NECZ
  284.  
  285. C     Passage à la partie dédiée selon la dimension courante
  286.       IF(IDIM.EQ.2) GO TO 20
  287.       IF(IDIM.EQ.3) GO TO 30
  288.  
  289. C     Sinon on provoque une erreur d'incompatibilité
  290. C     "Fonction indisponible en dimension %i1."
  291.       INTERR(1) = IDIM
  292.       CALL ERREUR(709)
  293.  
  294. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  295. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  296. CC
  297. CC  Dimension 2 (basé sur celui de dimension 3 avec Z en moins)
  298. CC
  299. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  300. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  301.  
  302. 20    CONTINUE
  303.  
  304. C     Recherche des bornes min/max des noeuds du maillage
  305. C     ***************************************************
  306.  
  307.       IND0 = (NUM(1,1)-1)*(IDIM+1)
  308.       XMIN0 = XCOOR(IND0+1)
  309.       YMIN0 = XCOOR(IND0+2)
  310.       XMAX0 = XCOOR(IND0+1)
  311.       YMAX0 = XCOOR(IND0+2)
  312.       DO 21 I = 1, NBN0
  313.           IND0 = (NUM(1,I)-1)*(IDIM+1)
  314.           X1 = XCOOR(IND0+1)
  315.           Y1 = XCOOR(IND0+2)
  316.           IF(X1.LT.XMIN0) XMIN0 = X1
  317.           IF(X1.GT.XMAX0) XMAX0 = X1
  318.           IF(Y1.LT.YMIN0) YMIN0 = Y1
  319.           IF(Y1.GT.YMAX0) YMAX0 = Y1
  320. 21    CONTINUE
  321.  
  322. C     Nombre de cellules suivant X et Y de la grille de cellules
  323. C     **********************************************************
  324.  
  325.       NX0 = INT((XMAX0-XMIN0)/C0) + 1
  326.       NY0 = INT((YMAX0-YMIN0)/C0) + 1
  327.  
  328. C     Répartition des noeuds dans les cellules (GRID2)
  329. C     ************************************************
  330.  
  331. C     Initialisation du segment contenant notre grille de cellules GRID2
  332. C     avec les dimensions NX0 et NY0 calculées juste avant
  333.       SEGINI GRID2SEG
  334. C     Boucle sur l'ensemble des noeuds
  335.       DO 22 I = 1, NBN0
  336. C         Récupération des coordonnées du noeud actuel
  337.           IND0 = (NUM(1,I)-1)*(IDIM+1)
  338.           X1 = XCOOR(IND0+1)
  339.           Y1 = XCOOR(IND0+2)
  340. C         Indices de la cellule contenant le noeud actuel
  341.           IX1 = INT((X1-XMIN0)/C0) + 1
  342.           IY1 = INT((Y1-YMIN0)/C0) + 1
  343. C         Ajout du noeud à la cellule correspondante
  344. C         Si la cellule n'existe pas encore
  345.           IF(GRID2(IX1,IY1).EQ.0) THEN
  346. C             Initialisation d'une liste d'entiers de taille 1
  347. C             des noeuds de cette cellule
  348.               JG = 1
  349.               SEGINI MLENTC
  350. C             On y ajoute le noeud actuel
  351.               MLENTC.LECT(1) = I
  352. C             On sauvegarde le pointeur de cette liste dans la grille
  353.               GRID2(IX1,IY1) = MLENTC
  354. C             On incrémente le nombre de cellule non vide
  355.               NBNEC0 = NBNEC0 + 1
  356. C             On ajoute les indices X et Y de la nouvelle cellule
  357. C             à la liste des indices des cellules non vides
  358.               JG = NBNEC0
  359.               SEGADJ NECX, NECY
  360.               NECX.LECT(JG) = IX1
  361.               NECY.LECT(JG) = IY1
  362.           ELSE
  363. C             On récupère le pointeur de la liste des noeuds de la cellule
  364.               MLENTC = GRID2(IX1,IY1)
  365. C             On ajuste la taille de cette liste
  366.               JG = MLENTC.LECT(/1) + 1
  367.               SEGADJ MLENTC
  368. C             On y ajoute le noeud actuel
  369.               MLENTC.LECT(JG) = I
  370.           ENDIF
  371. 22    CONTINUE
  372.  
  373. C     Calcul des indices X maximum des cellules distantes de R0 d'une
  374. C     cellule de référence d'indices X=0, Y=IMAX0+1 (INDEX2)
  375. C     ***************************************************************
  376.  
  377. C     Nombre maximum de cellules sur une distance de R0
  378.       IMAX0 = INT(R0/C0)
  379.       IF((R0/C0).GT.IMAX0) IMAX0 = IMAX0 + 1
  380. C     Cette même valeur + 1 (car souvent utilisée)
  381.       IMAX1 = IMAX0 + 1
  382. C     Initialisation du segment contenant les indices
  383. C     des cellules voisines à la cellule de référence
  384.       JG = 2*IMAX0 + 1
  385.       SEGINI INDEX2
  386.  
  387. C     La cellule la plus éloignée de la cellule de référence suivant X
  388. C     et dans un rayon R0 a pour coordonnées (IMAX0, IMAX0+1)
  389. C     On sauve donc IMAX0 à l'indice IMAX0+1 dans INDEX2
  390.       INDEX2.LECT(IMAX1) = IMAX0
  391.  
  392.  
  393. C     On parcourt les autres indices Y compris dans R0
  394. C     Boucle sur Y
  395.       DO 50 J = 1, IMAX0
  396. C         Calcul de l'indice X (IM1) de la cellule ayant pour indice
  397. C         Y=IMAX0+1+J la plus éloignée de la cellule de référence
  398. C         dans un rayon de R0
  399. C         x² + y² = r² => x = sqrt(r² - y²)
  400.           A0 = SQRT(R0**2 - ((J-1)*C0)**2)/C0
  401.           IM1 = INT(A0)
  402.           IF(A0.GT.IM1) IM1 = IM1 + 1
  403. C         Ajout de cette valeur à INDEX2 pour Y=IMAX0+1+J
  404.           INDEX2.LECT(IMAX1+J) = IM1
  405. C         La valeur de IM1 reste la même pour Y=IMAX0+1-J (symétrie)
  406.           INDEX2.LECT(IMAX1-J) = IM1
  407. 50    CONTINUE
  408.  
  409. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  410. C     Calcul des poids du voisinage de chaque noeud
  411. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  412.  
  413. C     On parcourt la liste des cellules non vides
  414. C     *******************************************
  415.  
  416.       DO 23 I = 1, NBNEC0
  417.  
  418. C         Coordonnées de la cellule actuelle
  419.           IX0 = NECX.LECT(I)
  420.           IY0 = NECY.LECT(I)
  421.  
  422. C         On génère la liste des noeuds voisins (MLENTV)
  423. C         °°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°
  424.  
  425. C         Initialisation liste de voisins potentiels pour toute la cellule
  426.           SEGINI MLENTV
  427. C         On initialise le nombre de noeuds voisins
  428.           JG = 0
  429. C         Boucles sur l'indices Y des cellules voisines
  430.           DO 51 J = 1, 2*IMAX0+1
  431.               IY1 = IY0 + IMAX1 - J
  432. C             Si l'indice Y de cellule voisine est un indice valide
  433.               IF((IY1.GE.1).AND.(IY1.LE.NY0)) THEN
  434.                   IM1 = INDEX2.LECT(J)
  435. C                 Boucle l'indice X des cellules voisines
  436.                   DO 52 K = -IM1, IM1
  437.                       IX1 = IX0 + K
  438. C                     Si l'indice X de cellule voisine est un indice valide
  439.                       IF((IX1.GE.1).AND.(IX1.LE.NX0)) THEN
  440. C                         Pointeur de la liste des noeuds de la cellule voisine
  441.                           MLENTC = GRID2(IX1,IY1)
  442. C                         Si la cellule voisine contient des noeuds
  443.                           IF(MLENTC.NE.0) THEN
  444.                               NBPC1 = MLENTC.LECT(/1)
  445. C                             On ajoute ces noeuds aux noeuds voisins de la cellule actuelle
  446.                               NBPV0 = JG
  447.                               JG = JG + NBPC1
  448.                               SEGADJ MLENTV
  449.                               DO 53 L = 1, NBPC1
  450.                                   MLENTV.LECT(NBPV0+L) = MLENTC.LECT(L)
  451. 53                            CONTINUE
  452.                           ENDIF
  453.                       ENDIF
  454. 52                CONTINUE
  455.               ENDIF
  456. 51        CONTINUE
  457.  
  458. C         On parcourt les noeuds de la cellule actuelle
  459. C         °°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°
  460.  
  461. C         Pointeur de la liste des noeuds de la cellule actuelle
  462.           MLENTC = GRID2(IX0,IY0)
  463. C         Nombre de noeuds contenus dans la cellule actuelle
  464.           NBPC0 = MLENTC.LECT(/1)
  465. C         Boucle sur chaque noeud de la cellule actuelle
  466.           DO 54, J = 1, NBPC0
  467. C             Récupération du noeud actuel
  468.               INDM0 = MLENTC.LECT(J)
  469.               NO0 = NUM(1,INDM0)
  470.               IND0 = (NO0-1)*(IDIM+1)
  471. C             Récupération des coordonnées du noeud actuel
  472.               X0 = XCOOR(IND0+1)
  473.               Y0 = XCOOR(IND0+2)
  474. C             Coordonnées -/+ R0
  475.               XLD0 = X0-R0
  476.               XLU0 = X0+R0
  477.               YLD0 = Y0-R0
  478.               YLU0 = Y0+R0
  479.  
  480. C             On génère la liste des noeuds voisins du noeud (MLENTN)
  481. C             du noeud actuel ainsi que leurs poids associés (MLREEP)
  482. C             -------------------------------------------------------
  483.  
  484. C             Initialisation de la liste des noeuds voisins (MLENTN)
  485. C             et de la liste de leurs poids respectifs (MLREEP)
  486.               JG = 1
  487.               SEGINI MLENTN, MLREEP
  488. C             On ajoute le noeud lui même et son poids propre (pour qu'il soit en premier dans la liste d'éléments pour que le dual corresponde)
  489.               MLENTN.LECT(1) = NO0
  490.               SELFW0 = VPOCHA(INDM0,1)
  491.               MLREEP.PROG(1) = SELFW0
  492. C             Initialisation de la somme des poids des noeuds voisins
  493.               SUMFAC0 = SELFW0
  494. C             On parcourt la liste des voisins potentiels
  495.               DO 24 K = 1, MLENTV.LECT(/1)
  496. C                 Récupération du noeud voisin
  497.                   INDM1 = MLENTV.LECT(K)
  498. C                 Si le noeud voisin traité est différent du noeud actuel
  499.                   IF(INDM1.NE.INDM0) THEN
  500.                       NO1 = NUM(1,INDM1)
  501.                       IND1 = (NO1-1)*(IDIM+1)
  502. C                     Coordonnées X du noeud voisin
  503.                       X1 = XCOOR(IND1+1)
  504. C                     Si la différence de coordonnées n'excède pas R0
  505.                       IF(X1.GE.XLD0.AND.X1.LE.XLU0) THEN
  506. C                         Coordonnées Y du noeud voisin
  507.                           Y1 = XCOOR(IND1+2)
  508. C                         Si la différence n'excède pas R0
  509.                           IF(Y1.GE.YLD0.AND.Y1.LE.YLU0) THEN
  510.                               A0 = R0 - SQRT((X0-X1)**2 +
  511.      &                                       (Y0-Y1)**2)
  512. C                             Si la distance n'excède pas R0
  513.                               IF(A0.GT.0.0) THEN
  514. C                                 Calcul du poids pour ce voisin
  515.                                   FAC0 = (A0/R0)**EMU0
  516.      &                                   * VPOCHA(INDM1,1)
  517. C                                 On sauvegarde le noeud voisin
  518. C                                 à la liste des noeuds voisins
  519. C                                 ainsi que son poids associé
  520.                                   JG = JG + 1
  521.                                   SEGADJ MLENTN, MLREEP
  522.                                   MLENTN.LECT(JG) = NO1
  523.                                   MLREEP.PROG(JG) = FAC0
  524. C                                 On met à jour la somme des poids
  525.                                   SUMFAC0 = SUMFAC0 + FAC0
  526.                               ENDIF
  527.                           ENDIF
  528.                       ENDIF
  529.                   ENDIF
  530. 24            CONTINUE
  531.  
  532. C             On ne conserve que les noeuds voisins de poids supérieur au critère
  533. C             -------------------------------------------------------------------
  534.  
  535. C             On boucle sur les noeuds voisins
  536.               A0 = CRI0 * SUMFAC0
  537.               DO 26 K = JG, 2, -1
  538.                   IF(MLREEP.PROG(K).LE.A0) THEN
  539. C                     On retire de la somme des poids
  540. C                     le poids du noeud voisin à éliminer
  541.                       SUMFAC0 = SUMFAC0 - MLREEP.PROG(K)
  542. C                     On élimine le noeud voisin et le poids associé
  543. C                      IF(K.NE.JG) THEN >> ce test est inutile car la boucle est ignorée si K>JG
  544.                           DO 27 L = K, JG-1
  545.                               MLENTN.LECT(L) = MLENTN.LECT(L+1)
  546.                               MLREEP.PROG(L) = MLREEP.PROG(L+1)
  547. 27                        CONTINUE
  548. C                      ENDIF
  549. C                     On met à jour le nombre de noeuds voisins
  550.                       JG = JG - 1
  551. C                     On ajuste la taille des listes
  552.                       SEGADJ MLENTN, MLREEP
  553.                   ENDIF
  554. 26            CONTINUE
  555.  
  556. C             Ajout d'une matrice élémentaire à la matrice de rigidité
  557. C             --------------------------------------------------------
  558.  
  559. C             On définit le nombre d'inconnues primales (= nombre de voisins)
  560.               NLIGRP = JG
  561. C             On définit le nombre de noeuds de l'élément de la matrice élémentaire
  562.               NBNN = JG
  563.  
  564.               IF(JG.GT.NBVMAX0) THEN
  565. C                 Si besoin on ajuste la taille de MLENTS
  566.                   NBVMAX0 = JG
  567.                   SEGADJ,MLENTS
  568.               ENDIF
  569.  
  570. C             Les matrices élémentaires de même nombre de noeuds (ou voisins) sont regroupées dans une même sous-matrice.
  571. C             Si la sous-matrice pour ce nombre de voisins n'existe pas encore, on l'initialise
  572.               IF((JG.GT.NBVMAX0).OR.(MLENTS.LECT(JG).EQ.0)) THEN
  573. C                 On met à jour le nombre de sous matrices
  574.                   NRIGEL = NRIGEL + 1
  575. C                 Sauvegarde de l'indice de cette nouvelle sous-matrice
  576.                   MLENTS.LECT(JG) = NRIGEL
  577. C                 Initialisation du segment descripteur de la sous-matrice
  578. C                 On indique son nombre d'inconnues duales (déjà fait au début)
  579. C                  NLIGRD = 1
  580.                   SEGINI DESCR
  581. C                 La première et unique inconnue duale correspondant au 1er noeud voisin (le noeud actuel)
  582.                   NOELED(1) = 1
  583. C                 Nom de l'inconnue dual
  584.                   LISDUA(1) = NOMD0
  585. C                 Les inconnues primales, allant simplement de 1 au nombre de voisins,
  586. C                 portent toutes le même nom
  587.                   DO 28 L = 1, NLIGRP
  588.                       NOELEP(L) = L
  589.                       LISINC(L) = NOMP0
  590. 28                CONTINUE
  591. C                 Initialisation d'un maillage
  592.                   NBELEM = 1
  593.                   NBSOUS = 0
  594.                   NBREF = 0
  595.                   SEGINI IPT1
  596. C                 Initialisation de la sous-matrice
  597. C                 Elle contiendra au moins 1 élément, qui va être ajouté
  598.                   NELRIG = 1
  599.                   SEGINI XMATRI
  600. C                 On indique que la matrice ne possède pas de symétries
  601.                   SYMRE = 2
  602. C                 On ajuste le nombre de sous matrices
  603.                   JG2 = JG
  604.                   JG = NRIGEL
  605.                   SEGADJ MLENTD, MLENTM, MLENTX
  606.                   JG = JG2
  607. C                 Sauvegarde des pointeurs des 3 segments crées ci-dessus
  608.                   MLENTD.LECT(NRIGEL) = DESCR
  609.                   MLENTM.LECT(NRIGEL) = IPT1
  610.                   MLENTX.LECT(NRIGEL) = XMATRI
  611. C             Si la sous-matrice existe, on la récupère
  612.               ELSE
  613.                   IPT1   = MLENTM.LECT(MLENTS.LECT(JG))
  614.                   XMATRI = MLENTX.LECT(MLENTS.LECT(JG))
  615. C                 On incrémente le nombre d'éléments du maillage
  616. C                 et le nombre de matrices élémentaires de la sous-matrice
  617.                   NBELEM = IPT1.NUM(/2) + 1
  618.                   NELRIG = RE(/3) + 1
  619.                   SEGADJ IPT1, XMATRI
  620.               ENDIF
  621. C             On boucle sur les noeuds voisins
  622.               DO 29 K = 1, JG
  623. C                 Dans IPT1, on ajoute la liste de noeuds voisins au maillage
  624.                   IPT1.NUM(K,NBELEM) = MLENTN.LECT(K)
  625. C                 Dans XMATRI, on ajoute la liste des composantes de la nouvelle matrice de rigidité élémentaire
  626.                   RE(1,K,NELRIG) = MLREEP.PROG(K)
  627. 29            CONTINUE
  628.  
  629. 54        CONTINUE
  630. 23    CONTINUE
  631. C     On supprime les segments devenus inutiles
  632.       SEGSUP GRID2SEG, INDEX2
  633.  
  634.       GO TO 98
  635.  
  636. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  637. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  638. CC
  639. CC  Dimension 3
  640. CC
  641. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  642. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  643.  
  644. 30    CONTINUE
  645.  
  646. C     Recherche des bornes min/max des noeuds du maillage
  647. C     ***************************************************
  648.  
  649.       IND0 = (NUM(1,1)-1)*(IDIM+1)
  650.       XMIN0 = XCOOR(IND0+1)
  651.       YMIN0 = XCOOR(IND0+2)
  652.       ZMIN0 = XCOOR(IND0+3)
  653.       XMAX0 = XCOOR(IND0+1)
  654.       YMAX0 = XCOOR(IND0+2)
  655.       ZMAX0 = XCOOR(IND0+3)
  656.       DO 31 I = 2, NBN0
  657.           IND0 = (NUM(1,I)-1)*(IDIM+1)
  658.           X1 = XCOOR(IND0+1)
  659.           Y1 = XCOOR(IND0+2)
  660.           Z1 = XCOOR(IND0+3)
  661.           IF(X1.LT.XMIN0) XMIN0 = X1
  662.           IF(X1.GT.XMAX0) XMAX0 = X1
  663.           IF(Y1.LT.YMIN0) YMIN0 = Y1
  664.           IF(Y1.GT.YMAX0) YMAX0 = Y1
  665.           IF(Z1.LT.ZMIN0) ZMIN0 = Z1
  666.           IF(Z1.GT.ZMAX0) ZMAX0 = Z1
  667. 31    CONTINUE
  668.  
  669. C     Nombre de cellules suivant X, Y et Z de la grille de cellules
  670. C     *************************************************************
  671.  
  672.       NX0 = INT((XMAX0-XMIN0)/C0) + 1
  673.       NY0 = INT((YMAX0-YMIN0)/C0) + 1
  674.       NZ0 = INT((ZMAX0-ZMIN0)/C0) + 1
  675.  
  676. C     Répartition des noeuds dans les cellules (GRID2)
  677. C     ************************************************
  678.  
  679. C     Initialisation du segment contenant notre grille de cellules GRID3
  680. C     avec les dimensions NX0, NY0 et NZ0 calculées juste avant
  681.       SEGINI GRID3SEG
  682. C     Boucle sur l'ensemble des noeuds
  683.       DO 32 I = 1, NBN0
  684. C         Récupération des coordonnées du noeud actuel
  685.           IND0 = (NUM(1,I)-1)*(IDIM+1)
  686.           X1 = XCOOR(IND0+1)
  687.           Y1 = XCOOR(IND0+2)
  688.           Z1 = XCOOR(IND0+3)
  689. C         Indices de la cellule contenant le noeud actuel
  690.           IX1 = INT((X1-XMIN0)/C0) + 1
  691.           IY1 = INT((Y1-YMIN0)/C0) + 1
  692.           IZ1 = INT((Z1-ZMIN0)/C0) + 1
  693. C         Ajout du noeud à la cellule correspondante
  694. C         Si la cellule n'existe pas encore
  695.           IF(GRID3(IX1,IY1,IZ1).EQ.0) THEN
  696. C             Initialisation d'une liste d'entiers de taille 1
  697. C             des noeuds de cette cellule
  698.               JG = 1
  699.               SEGINI MLENTC
  700. C             On y ajoute le noeud actuel
  701.               MLENTC.LECT(1) = I
  702. C             On sauvegarde le pointeur de cette liste dans la grille
  703.               GRID3(IX1,IY1,IZ1) = MLENTC
  704. C             On incrémente le nombre de cellule non vide
  705.               NBNEC0 = NBNEC0 + 1
  706. C             On ajoute les indices X, Y et Z de la nouvelle cellule
  707. C             à la liste des indices des cellules non vides
  708.               JG = NBNEC0
  709.               SEGADJ NECX, NECY, NECZ
  710.               NECX.LECT(JG) = IX1
  711.               NECY.LECT(JG) = IY1
  712.               NECZ.LECT(JG) = IZ1
  713.           ELSE
  714. C             On récupère le pointeur de la liste des noeuds de la cellule
  715.               MLENTC = GRID3(IX1,IY1,IZ1)
  716. C             On ajuste la taille de cette liste
  717.               JG = MLENTC.LECT(/1) + 1
  718.               SEGADJ MLENTC
  719. C             On y ajoute le noeud actuel
  720.               MLENTC.LECT(JG) = I
  721.           ENDIF
  722. 32    CONTINUE
  723.  
  724. C     Calcul des indices X maximum des cellules distantes de R0 d'une
  725. C     cellule de référence d'indices X=0, Y=IMAX0+1, Z=IMAX0+1 (INDEX3)
  726. C     *****************************************************************
  727.  
  728. C     Nombre maximum de cellules sur une distance de R0
  729.       IMAX0 = INT(R0/C0)
  730.       IF((R0/C0).GT.IMAX0) IMAX0 = IMAX0 + 1
  731. C     Cette même valeur + 1 (car souvent utilisée)
  732.       IMAX1 = IMAX0 + 1
  733. C     Initialisation du segment contenant les indices
  734. C     des cellules voisines à la cellule de référence
  735.       JG = 2*IMAX0 + 1
  736.       SEGINI INDEX3SEG
  737.  
  738. C     La cellule la plus éloignée de la cellule de référence suivant X
  739. C     et dans un rayon R0 a pour coordonnées (IMAX0, IMAX0+1, IMAX0+1)
  740. C     On sauve donc IMAX0 à l'indice IMAX0+1, IMAX0+1 dans INDEX3
  741.       INDEX3(IMAX1,IMAX1) = IMAX0
  742.  
  743. C     On parcourt les autres indices Y et Z compris dans R0
  744. C     Boucle sur Y
  745.       DO 60 J = 1, IMAX0
  746. C         Calcul de l'indice X (IM1) de la cellule ayant pour indice
  747. C         Y=IMAX0+1+J et pour indice Z=IMAX0+1 la plus éloignée de la
  748. C         cellule de référence dans un rayon de R0
  749. C         x² + y² = r² => x = sqrt(r² - y²), pour z = 0
  750.           A0 = SQRT(R0**2 - ((J-1)*C0)**2)/C0
  751.           IM1 = INT(A0)
  752.           IF(A0.GT.IM1) IM1 = IM1 + 1
  753. C         Ajout de cette valeur à INDEX3 pour Y=IMAX0+1+J et Z=IMAX0+1
  754.           INDEX3(IMAX1+J,IMAX1) = IM1
  755. C         La valeur de IM1 reste la même pour Y=IMAX0+1-J (symétrie)
  756.           INDEX3(IMAX1-J,IMAX1) = IM1
  757. C         Enfin la valeur de IM1 reste la même quand on inverse Y et Z
  758.           INDEX3(IMAX1,IMAX1+J) = IM1
  759.           INDEX3(IMAX1,IMAX1-J) = IM1
  760. C         Boucle sur Z
  761.           DO 61 K = 1, IM1
  762. C             Calcul de l'indice X (IM2) de la cellule ayant pour
  763. C             indice Y=IMAX0+1+J et pour indice Z=IMAX0+1+K la plus
  764. C             éloignée de la cellule de référence dans un rayon de R0
  765. C             x² + y² + z² = r² => x = sqrt(r² - y² - z²)
  766.               A0 = SQRT(R0**2 - ((J-1)*C0)**2 - ((K-1)*C0)**2)/C0
  767.               IM2 = INT(A0)
  768.               IF(A0.GT.IM2) IM2 = IM2 + 1
  769. C             Ajout du cette valeur à INDEX3
  770.               INDEX3(IMAX1+J,IMAX1+K) = IM2
  771. C             Symétrie
  772.               INDEX3(IMAX1+J,IMAX1-K) = IM2
  773. C             Inversion de Y et Z
  774.               INDEX3(IMAX1-J,IMAX1+K) = IM2
  775.               INDEX3(IMAX1-J,IMAX1-K) = IM2
  776. 61        CONTINUE
  777. 60    CONTINUE
  778.  
  779. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  780. C     Calcul des poids du voisinage de chaque noeud
  781. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  782.  
  783. C     On parcourt la liste des cellules non vides
  784. C     *******************************************
  785.  
  786.       DO 33 I = 1, NBNEC0
  787.  
  788. C         Indices de la cellule actuelle
  789.           IX0 = NECX.LECT(I)
  790.           IY0 = NECY.LECT(I)
  791.           IZ0 = NECZ.LECT(I)
  792.  
  793. C         On génère la liste des noeuds voisins (MLENTV)
  794. C         °°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°
  795.  
  796. C         Initialisation liste de voisins potentiels de la cellule actuelle
  797.           SEGINI MLENTV
  798. C         On initialise le nombre de noeuds voisins
  799.           JG = 0
  800. C         Boucles sur l'indices Y des cellules voisines
  801.           DO 62 J = 1, 2*IMAX0+1
  802.               IY1 = IY0 + IMAX1 - J
  803. C             Si l'indice Y de cellule voisine est un indice valide
  804.               IF((IY1.GE.1).AND.(IY1.LE.NY0)) THEN
  805.                   IM1 = INDEX3(J,IMAX1)
  806. C                 Boucle l'indice Z des cellules voisines
  807.                   DO 63 K = IMAX1-IM1, IMAX1+IM1
  808.                       IZ1 = IZ0 + IMAX1 - K
  809. C                     Si l'indice Z de cellule voisine est un indice valide
  810.                       IF((IZ1.GE.1).AND.(IZ1.LE.NZ0)) THEN
  811.                           IM2 = INDEX3(J,K)
  812. C                         Boucle l'indice X des cellules voisines
  813.                           DO 64 L = -IM2, IM2
  814.                               IX1 = IX0 + L
  815. C                             Si l'indice X de cellule voisine est un indice valide
  816.                               IF((IX1.GE.1).AND.(IX1.LE.NX0)) THEN
  817. C                                 Pointeur de la liste des noeuds de la cellule voisine
  818.                                   MLENTC = GRID3(IX1,IY1,IZ1)
  819. C                                 Si la cellule voisine contient des noeuds
  820.                                   IF(MLENTC.NE.0) THEN
  821.                                       NBPC0 = MLENTC.LECT(/1)
  822. C                                     On ajoute ces noeuds aux noeuds voisins de la cellule actuelle
  823.                                       NBPV0 = JG
  824.                                       JG = JG + NBPC0
  825.                                       SEGADJ MLENTV
  826.                                       DO 65 M = 1, NBPC0
  827.                                           MLENTV.LECT(NBPV0+M) =
  828.      &                                        MLENTC.LECT(M)
  829. 65                                    CONTINUE
  830.                                   ENDIF
  831.                               ENDIF
  832. 64                        CONTINUE
  833.                       ENDIF
  834. 63                CONTINUE
  835.               ENDIF
  836. 62        CONTINUE
  837.  
  838. C         On parcourt les noeuds de la cellule actuelle
  839. C         °°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°
  840.  
  841. C         Pointeur de la liste des noeuds de la cellule actuelle
  842.           MLENTC = GRID3(IX0,IY0,IZ0)
  843. C         Nombre de noeuds contenus dans la cellule actuelle
  844.           NBPC0 = MLENTC.LECT(/1)
  845. C         Boucle sur chaque noeud de la cellule actuelle
  846.           DO 66, J = 1, NBPC0
  847. C             Récupération du noeud actuel
  848.               INDM0 = MLENTC.LECT(J)
  849.               NO0 = NUM(1,INDM0)
  850.               IND0 = (NO0-1)*(IDIM+1)
  851. C             Récupération des coordonnées du noeud actuel
  852.               X0 = XCOOR(IND0+1)
  853.               Y0 = XCOOR(IND0+2)
  854.               Z0 = XCOOR(IND0+3)
  855. C             Coordonnées -/+ R0
  856.               XLD0 = X0-R0
  857.               XLU0 = X0+R0
  858.               YLD0 = Y0-R0
  859.               YLU0 = Y0+R0
  860.               ZLD0 = Z0-R0
  861.               ZLU0 = Z0+R0
  862.  
  863. C             On génère la liste des noeuds voisins du noeud (MLENTN)
  864. C             du noeud actuel ainsi que leurs poids associés (MLREEP)
  865. C             -------------------------------------------------------
  866.  
  867. C             Initialisation de la liste des noeuds voisins (MLENTN)
  868. C             et de la liste de leurs poids respectifs (MLREEP)
  869.               JG = 1
  870.               SEGINI MLENTN, MLREEP
  871. C             On ajoute le noeud lui même et son  poids propre (pour qu'il soit en premier dans la liste d'éléments pour que le dual corresponde)
  872.               MLENTN.LECT(1) = NO0
  873.               SELFW0 = VPOCHA(INDM0,1)
  874.               MLREEP.PROG(1) = SELFW0
  875. C             Initialisation de la somme des poids des noeuds voisins
  876.               SUMFAC0 = SELFW0
  877. C             On parcourt la liste des voisins potentiels
  878.               DO 34 K = 1, MLENTV.LECT(/1)
  879. C                 Récupération du noeud voisin
  880.                   INDM1 = MLENTV.LECT(K)
  881. C                 Si le noeud voisin traité est différent du noeud actuel
  882.                   IF(INDM1.NE.INDM0) THEN
  883.                       NO1 = NUM(1,INDM1)
  884.                       IND1 = (NO1-1)*(IDIM+1)
  885. C                     Coordonnées X du noeud voisin
  886.                       X1 = XCOOR(IND1+1)
  887. C                     Si la différence de coordonnées n'excède pas R0
  888.                       IF(X1.GE.XLD0.AND.X1.LE.XLU0) THEN
  889. C                         Coordonnées Y du noeud voisin
  890.                           Y1 = XCOOR(IND1+2)
  891. C                         Si la différence n'excède pas R0
  892.                           IF(Y1.GE.YLD0.AND.Y1.LE.YLU0) THEN
  893. C                             Coordonnées Z du noeud voisin
  894.                               Z1 = XCOOR(IND1+3)
  895. C                             Si la différence n'excède pas R0
  896.                               IF(Z1.GE.ZLD0.AND.Z1.LE.ZLU0) THEN
  897.                                   A0 = R0 - SQRT((X0-X1)**2 +
  898.      &                                       (Y0-Y1)**2 + (Z0-Z1)**2)
  899. C                                 Si la distance n'excède pas R0
  900.                                   IF(A0.GT.0.0) THEN
  901. C                                     Calcul du poids pour ce voisin
  902.                                       FAC0 = (A0/R0)**EMU0
  903.      &                                       * VPOCHA(INDM1,1)
  904. C                                     On sauvegarde le noeud voisin
  905. C                                     à la liste des noeuds voisins
  906. C                                     ainsi que son poids associé
  907.                                       JG = JG + 1
  908.                                       SEGADJ MLENTN, MLREEP
  909.                                       MLENTN.LECT(JG) = NO1
  910.                                       MLREEP.PROG(JG) = FAC0
  911. C                                     On met à jour la somme des poids
  912.                                       SUMFAC0 = SUMFAC0 + FAC0
  913.                                   ENDIF
  914.                               ENDIF
  915.                           ENDIF
  916.                       ENDIF
  917.                   ENDIF
  918. 34            CONTINUE
  919.  
  920. C             On ne conserve que les noeuds voisins de poids supérieur au critère
  921. C             -------------------------------------------------------------------
  922.  
  923. C             On boucle sur les noeuds voisins
  924.               A0 = CRI0 * SUMFAC0
  925.               DO 36 K = JG, 2, -1
  926.                   IF(MLREEP.PROG(K).LE.A0) THEN
  927. C                     On retire de la somme des poids
  928. C                     le poids du noeud voisin à éliminer
  929.                       SUMFAC0 = SUMFAC0 - MLREEP.PROG(K)
  930. C                     On élimine le noeud voisin et le poids associé
  931. C                      IF(K.NE.JG) THEN >> ce test est inutile car la boucle est ignorée si K>JG
  932.                           DO 38 L = K, JG-1
  933.                               MLENTN.LECT(L) = MLENTN.LECT(L+1)
  934.                               MLREEP.PROG(L) = MLREEP.PROG(L+1)
  935. 38                        CONTINUE
  936. C                      ENDIF
  937. C                     On met à jour le nombre de noeuds voisins
  938.                       JG = JG - 1
  939. C                     On ajuste la taille des listes
  940.                       SEGADJ MLENTN, MLREEP
  941.                   ENDIF
  942. 36            CONTINUE
  943.  
  944. C             Ajout d'une matrice élémentaire à la matrice de rigidité
  945. C             --------------------------------------------------------
  946.  
  947. C             On définit le nombre d'inconnues primales (= nombre de voisins)
  948.               NLIGRP = JG
  949. C             On définit le nombre de noeuds de l'élément de la matrice élémentaire
  950.               NBNN = JG
  951.  
  952.               IF(JG.GT.NBVMAX0) THEN
  953. C                 Si besoin on ajuste la taille de MLENTS
  954.                   NBVMAX0 = JG
  955.                   SEGADJ MLENTS
  956.               ENDIF
  957.  
  958. C             Les matrices élémentaires de même nombre de noeuds (ou voisins) sont regroupées dans une même sous-matrice.
  959. C             Si la sous-matrice pour ce nombre de voisins n'existe pas encore, on l'initialise
  960.               IF((JG.GT.NBVMAX0).OR.(MLENTS.LECT(JG).EQ.0)) THEN
  961. C                 On met à jour le nombre de sous matrices
  962.                   NRIGEL = NRIGEL + 1
  963. C                 Sauvegarde de l'indice de cette nouvelle sous-matrice
  964.                   MLENTS.LECT(JG) = NRIGEL
  965. C                 Initialisation du segment descripteur de la sous-matrice
  966. C                 On indique son nombre d'inconnues duales (déjà fait au début)
  967. C                  NLIGRD = 1
  968.                   SEGINI DESCR
  969. C                 La première et unique inconnue duale correspondant au 1er noeud voisin (le noeud actuel)
  970.                   NOELED(1) = 1
  971. C                 Nom de l'inconnue dual
  972.                   LISDUA(1) = NOMD0
  973. C                 Les inconnues primales, allant simplement de 1 au nombre de voisins,
  974. C                 portent toutes le même nom
  975.                   DO 39 L = 1, NLIGRP
  976.                       NOELEP(L) = L
  977.                       LISINC(L) = NOMP0
  978. 39                CONTINUE
  979. C                 Initialisation d'un maillage
  980.                   NBELEM = 1
  981.                   NBSOUS = 0
  982.                   NBREF = 0
  983.                   SEGINI IPT1
  984. C                 Initialisation de la sous-matrice
  985. C                 Elle contiendra au moins 1 élément, qui va être ajouté
  986.                   NELRIG = 1
  987.                   SEGINI XMATRI
  988. C                 On indique que la matrice ne possède pas de symétries
  989.                   SYMRE = 2
  990. C                 On ajuste le nombre de sous matrices
  991.                   JG2 = JG
  992.                   JG = NRIGEL
  993.                   SEGADJ MLENTD, MLENTM, MLENTX
  994.                   JG = JG2
  995. C                 Sauvegarde des pointeurs des 3 segments crées ci-dessus
  996.                   MLENTD.LECT(NRIGEL) = DESCR
  997.                   MLENTM.LECT(NRIGEL) = IPT1
  998.                   MLENTX.LECT(NRIGEL) = XMATRI
  999. C             Si la sous-matrice existe, on la récupère
  1000.               ELSE
  1001.                   IPT1 = MLENTM.LECT(MLENTS.LECT(JG))
  1002.                   XMATRI = MLENTX.LECT(MLENTS.LECT(JG))
  1003. C                 On incrémente le nombre d'éléments du maillage
  1004. C                 et le nombre de matrices élémentaires de la sous-matrice
  1005.                   NBELEM = IPT1.NUM(/2) + 1
  1006.                   NELRIG = RE(/3) + 1
  1007.                   SEGADJ IPT1, XMATRI
  1008.               ENDIF
  1009. C             On boucle sur les noeuds voisins
  1010.               DO 40 K = 1, JG
  1011. C                 Dans IPT1, on ajoute la liste de noeuds voisins au maillage
  1012.                   IPT1.NUM(K,NBELEM) = MLENTN.LECT(K)
  1013. C                 Dans XMATRI, on ajoute la liste des composantes de la nouvelle matrice de rigidité élémentaire
  1014.                   RE(1,K,NELRIG) = MLREEP.PROG(K)
  1015. 40            CONTINUE
  1016. 66        CONTINUE
  1017. 33    CONTINUE
  1018. C     On supprime les segments devenus inutiles
  1019.       SEGSUP GRID3SEG, INDEX3SEG
  1020.  
  1021. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1022. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1023. CC
  1024. CC  Combinaison des sous-matrices dans une matrice globale de rigidité
  1025. CC
  1026. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1027. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1028.  
  1029. 98    CONTINUE
  1030.  
  1031. C     On initialise le segment de la matrice de rigidité globale
  1032. C     qui prendra la taille de NRIGEL
  1033.       SEGINI MRIGID
  1034. C     On indique son type de matrice
  1035.       MTYMAT = 'RIGIDITE'
  1036.       IFORIG = mchpoi.IFOPOI
  1037.  
  1038. C     On remplit la matrice de rigidité avec les matrices élémentaires,
  1039. C     descripteurs et maillages correspondant à chaque sous-matrice
  1040.       DO 11 I = 1, NRIGEL
  1041. C               Récupération des segments du descripteur,
  1042. C               de la sous-matrice et de la liste des noeuds
  1043.                 DESCR  = MLENTD.LECT(I)
  1044.                 XMATRI = MLENTX.LECT(I)
  1045.                 IPT1   = MLENTM.LECT(I)
  1046. C               Affectation du maillage, du descripteur et du contenu de la sous-matrice
  1047.                 MRIGID.IRIGEL(1,I) = IPT1
  1048.                 MRIGID.IRIGEL(3,I) = DESCR
  1049.                 MRIGID.IRIGEL(4,I) = XMATRI
  1050. C               Retrait du caractère *MOD des SEGMENTS crees dans ce SP                
  1051.                 SEGACT,IPT1,DESCR,XMATRI
  1052. C               On indique que la matrice globale ne possède pas de symétries
  1053.                 MRIGID.IRIGEL(7,I) = 2
  1054. C               On précise son coefficient multiplicateur
  1055.                 MRIGID.COERIG(I) = 1.D0
  1056. 11    CONTINUE
  1057.       SEGACT,MRIGID
  1058.  
  1059. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1060. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1061. CC
  1062. CC  Nettoyage final et sortie
  1063. CC
  1064. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1065. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1066.  
  1067. C     On supprime les segments utiles seulement au sein de la subroutine et quand R0 != 0
  1068.       SEGSUP MLENTX, MLENTD, MLENTM, MLENTN, MLREEP, NECX, NECY, NECZ
  1069.  
  1070. 99    CONTINUE
  1071.  
  1072.       SEGDES MCOORD
  1073.  
  1074. C     On désactive les segments restants avant de sortir
  1075.       SEGDES MRIGID, MELEME, MPOVAL, MSOUPO, MCHPOI
  1076. C     On renvoie la matrice de rigidité
  1077.       CALL ECROBJ('RIGIDITE',MRIGID)
  1078.  
  1079.       RETURN
  1080.       END
  1081.  
  1082.  

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