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    PV090527  26/04/28    21:15:55     12529          
  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.                   rigrel=0
  600.                   SEGINI XMATRI
  601. C                 On indique que la matrice ne possède pas de symétries
  602.                   SYMRE = 2
  603. C                 On ajuste le nombre de sous matrices
  604.                   JG2 = JG
  605.                   JG = NRIGEL
  606.                   SEGADJ MLENTD, MLENTM, MLENTX
  607.                   JG = JG2
  608. C                 Sauvegarde des pointeurs des 3 segments crées ci-dessus
  609.                   MLENTD.LECT(NRIGEL) = DESCR
  610.                   MLENTM.LECT(NRIGEL) = IPT1
  611.                   MLENTX.LECT(NRIGEL) = XMATRI
  612. C             Si la sous-matrice existe, on la récupère
  613.               ELSE
  614.                   IPT1   = MLENTM.LECT(MLENTS.LECT(JG))
  615.                   XMATRI = MLENTX.LECT(MLENTS.LECT(JG))
  616. C                 On incrémente le nombre d'éléments du maillage
  617. C                 et le nombre de matrices élémentaires de la sous-matrice
  618.                   NBELEM = IPT1.NUM(/2) + 1
  619.                   NELRIG = RE(/3) + 1
  620.                   rigrel=0
  621.                   SEGADJ IPT1, XMATRI
  622.               ENDIF
  623. C             On boucle sur les noeuds voisins
  624.               DO 29 K = 1, JG
  625. C                 Dans IPT1, on ajoute la liste de noeuds voisins au maillage
  626.                   IPT1.NUM(K,NBELEM) = MLENTN.LECT(K)
  627. C                 Dans XMATRI, on ajoute la liste des composantes de la nouvelle matrice de rigidité élémentaire
  628.                   RE(1,K,NELRIG) = MLREEP.PROG(K)
  629. 29            CONTINUE
  630.  
  631. 54        CONTINUE
  632. 23    CONTINUE
  633. C     On supprime les segments devenus inutiles
  634.       SEGSUP GRID2SEG, INDEX2
  635.  
  636.       GO TO 98
  637.  
  638. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  639. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  640. CC
  641. CC  Dimension 3
  642. CC
  643. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  644. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  645.  
  646. 30    CONTINUE
  647.  
  648. C     Recherche des bornes min/max des noeuds du maillage
  649. C     ***************************************************
  650.  
  651.       IND0 = (NUM(1,1)-1)*(IDIM+1)
  652.       XMIN0 = XCOOR(IND0+1)
  653.       YMIN0 = XCOOR(IND0+2)
  654.       ZMIN0 = XCOOR(IND0+3)
  655.       XMAX0 = XCOOR(IND0+1)
  656.       YMAX0 = XCOOR(IND0+2)
  657.       ZMAX0 = XCOOR(IND0+3)
  658.       DO 31 I = 2, NBN0
  659.           IND0 = (NUM(1,I)-1)*(IDIM+1)
  660.           X1 = XCOOR(IND0+1)
  661.           Y1 = XCOOR(IND0+2)
  662.           Z1 = XCOOR(IND0+3)
  663.           IF(X1.LT.XMIN0) XMIN0 = X1
  664.           IF(X1.GT.XMAX0) XMAX0 = X1
  665.           IF(Y1.LT.YMIN0) YMIN0 = Y1
  666.           IF(Y1.GT.YMAX0) YMAX0 = Y1
  667.           IF(Z1.LT.ZMIN0) ZMIN0 = Z1
  668.           IF(Z1.GT.ZMAX0) ZMAX0 = Z1
  669. 31    CONTINUE
  670.  
  671. C     Nombre de cellules suivant X, Y et Z de la grille de cellules
  672. C     *************************************************************
  673.  
  674.       NX0 = INT((XMAX0-XMIN0)/C0) + 1
  675.       NY0 = INT((YMAX0-YMIN0)/C0) + 1
  676.       NZ0 = INT((ZMAX0-ZMIN0)/C0) + 1
  677.  
  678. C     Répartition des noeuds dans les cellules (GRID2)
  679. C     ************************************************
  680.  
  681. C     Initialisation du segment contenant notre grille de cellules GRID3
  682. C     avec les dimensions NX0, NY0 et NZ0 calculées juste avant
  683.       SEGINI GRID3SEG
  684. C     Boucle sur l'ensemble des noeuds
  685.       DO 32 I = 1, NBN0
  686. C         Récupération des coordonnées du noeud actuel
  687.           IND0 = (NUM(1,I)-1)*(IDIM+1)
  688.           X1 = XCOOR(IND0+1)
  689.           Y1 = XCOOR(IND0+2)
  690.           Z1 = XCOOR(IND0+3)
  691. C         Indices de la cellule contenant le noeud actuel
  692.           IX1 = INT((X1-XMIN0)/C0) + 1
  693.           IY1 = INT((Y1-YMIN0)/C0) + 1
  694.           IZ1 = INT((Z1-ZMIN0)/C0) + 1
  695. C         Ajout du noeud à la cellule correspondante
  696. C         Si la cellule n'existe pas encore
  697.           IF(GRID3(IX1,IY1,IZ1).EQ.0) THEN
  698. C             Initialisation d'une liste d'entiers de taille 1
  699. C             des noeuds de cette cellule
  700.               JG = 1
  701.               SEGINI MLENTC
  702. C             On y ajoute le noeud actuel
  703.               MLENTC.LECT(1) = I
  704. C             On sauvegarde le pointeur de cette liste dans la grille
  705.               GRID3(IX1,IY1,IZ1) = MLENTC
  706. C             On incrémente le nombre de cellule non vide
  707.               NBNEC0 = NBNEC0 + 1
  708. C             On ajoute les indices X, Y et Z de la nouvelle cellule
  709. C             à la liste des indices des cellules non vides
  710.               JG = NBNEC0
  711.               SEGADJ NECX, NECY, NECZ
  712.               NECX.LECT(JG) = IX1
  713.               NECY.LECT(JG) = IY1
  714.               NECZ.LECT(JG) = IZ1
  715.           ELSE
  716. C             On récupère le pointeur de la liste des noeuds de la cellule
  717.               MLENTC = GRID3(IX1,IY1,IZ1)
  718. C             On ajuste la taille de cette liste
  719.               JG = MLENTC.LECT(/1) + 1
  720.               SEGADJ MLENTC
  721. C             On y ajoute le noeud actuel
  722.               MLENTC.LECT(JG) = I
  723.           ENDIF
  724. 32    CONTINUE
  725.  
  726. C     Calcul des indices X maximum des cellules distantes de R0 d'une
  727. C     cellule de référence d'indices X=0, Y=IMAX0+1, Z=IMAX0+1 (INDEX3)
  728. C     *****************************************************************
  729.  
  730. C     Nombre maximum de cellules sur une distance de R0
  731.       IMAX0 = INT(R0/C0)
  732.       IF((R0/C0).GT.IMAX0) IMAX0 = IMAX0 + 1
  733. C     Cette même valeur + 1 (car souvent utilisée)
  734.       IMAX1 = IMAX0 + 1
  735. C     Initialisation du segment contenant les indices
  736. C     des cellules voisines à la cellule de référence
  737.       JG = 2*IMAX0 + 1
  738.       SEGINI INDEX3SEG
  739.  
  740. C     La cellule la plus éloignée de la cellule de référence suivant X
  741. C     et dans un rayon R0 a pour coordonnées (IMAX0, IMAX0+1, IMAX0+1)
  742. C     On sauve donc IMAX0 à l'indice IMAX0+1, IMAX0+1 dans INDEX3
  743.       INDEX3(IMAX1,IMAX1) = IMAX0
  744.  
  745. C     On parcourt les autres indices Y et Z compris dans R0
  746. C     Boucle sur Y
  747.       DO 60 J = 1, IMAX0
  748. C         Calcul de l'indice X (IM1) de la cellule ayant pour indice
  749. C         Y=IMAX0+1+J et pour indice Z=IMAX0+1 la plus éloignée de la
  750. C         cellule de référence dans un rayon de R0
  751. C         x² + y² = r² => x = sqrt(r² - y²), pour z = 0
  752.           A0 = SQRT(R0**2 - ((J-1)*C0)**2)/C0
  753.           IM1 = INT(A0)
  754.           IF(A0.GT.IM1) IM1 = IM1 + 1
  755. C         Ajout de cette valeur à INDEX3 pour Y=IMAX0+1+J et Z=IMAX0+1
  756.           INDEX3(IMAX1+J,IMAX1) = IM1
  757. C         La valeur de IM1 reste la même pour Y=IMAX0+1-J (symétrie)
  758.           INDEX3(IMAX1-J,IMAX1) = IM1
  759. C         Enfin la valeur de IM1 reste la même quand on inverse Y et Z
  760.           INDEX3(IMAX1,IMAX1+J) = IM1
  761.           INDEX3(IMAX1,IMAX1-J) = IM1
  762. C         Boucle sur Z
  763.           DO 61 K = 1, IM1
  764. C             Calcul de l'indice X (IM2) de la cellule ayant pour
  765. C             indice Y=IMAX0+1+J et pour indice Z=IMAX0+1+K la plus
  766. C             éloignée de la cellule de référence dans un rayon de R0
  767. C             x² + y² + z² = r² => x = sqrt(r² - y² - z²)
  768.               A0 = SQRT(R0**2 - ((J-1)*C0)**2 - ((K-1)*C0)**2)/C0
  769.               IM2 = INT(A0)
  770.               IF(A0.GT.IM2) IM2 = IM2 + 1
  771. C             Ajout du cette valeur à INDEX3
  772.               INDEX3(IMAX1+J,IMAX1+K) = IM2
  773. C             Symétrie
  774.               INDEX3(IMAX1+J,IMAX1-K) = IM2
  775. C             Inversion de Y et Z
  776.               INDEX3(IMAX1-J,IMAX1+K) = IM2
  777.               INDEX3(IMAX1-J,IMAX1-K) = IM2
  778. 61        CONTINUE
  779. 60    CONTINUE
  780.  
  781. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  782. C     Calcul des poids du voisinage de chaque noeud
  783. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  784.  
  785. C     On parcourt la liste des cellules non vides
  786. C     *******************************************
  787.  
  788.       DO 33 I = 1, NBNEC0
  789.  
  790. C         Indices de la cellule actuelle
  791.           IX0 = NECX.LECT(I)
  792.           IY0 = NECY.LECT(I)
  793.           IZ0 = NECZ.LECT(I)
  794.  
  795. C         On génère la liste des noeuds voisins (MLENTV)
  796. C         °°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°
  797.  
  798. C         Initialisation liste de voisins potentiels de la cellule actuelle
  799.           SEGINI MLENTV
  800. C         On initialise le nombre de noeuds voisins
  801.           JG = 0
  802. C         Boucles sur l'indices Y des cellules voisines
  803.           DO 62 J = 1, 2*IMAX0+1
  804.               IY1 = IY0 + IMAX1 - J
  805. C             Si l'indice Y de cellule voisine est un indice valide
  806.               IF((IY1.GE.1).AND.(IY1.LE.NY0)) THEN
  807.                   IM1 = INDEX3(J,IMAX1)
  808. C                 Boucle l'indice Z des cellules voisines
  809.                   DO 63 K = IMAX1-IM1, IMAX1+IM1
  810.                       IZ1 = IZ0 + IMAX1 - K
  811. C                     Si l'indice Z de cellule voisine est un indice valide
  812.                       IF((IZ1.GE.1).AND.(IZ1.LE.NZ0)) THEN
  813.                           IM2 = INDEX3(J,K)
  814. C                         Boucle l'indice X des cellules voisines
  815.                           DO 64 L = -IM2, IM2
  816.                               IX1 = IX0 + L
  817. C                             Si l'indice X de cellule voisine est un indice valide
  818.                               IF((IX1.GE.1).AND.(IX1.LE.NX0)) THEN
  819. C                                 Pointeur de la liste des noeuds de la cellule voisine
  820.                                   MLENTC = GRID3(IX1,IY1,IZ1)
  821. C                                 Si la cellule voisine contient des noeuds
  822.                                   IF(MLENTC.NE.0) THEN
  823.                                       NBPC0 = MLENTC.LECT(/1)
  824. C                                     On ajoute ces noeuds aux noeuds voisins de la cellule actuelle
  825.                                       NBPV0 = JG
  826.                                       JG = JG + NBPC0
  827.                                       SEGADJ MLENTV
  828.                                       DO 65 M = 1, NBPC0
  829.                                           MLENTV.LECT(NBPV0+M) =
  830.      &                                        MLENTC.LECT(M)
  831. 65                                    CONTINUE
  832.                                   ENDIF
  833.                               ENDIF
  834. 64                        CONTINUE
  835.                       ENDIF
  836. 63                CONTINUE
  837.               ENDIF
  838. 62        CONTINUE
  839.  
  840. C         On parcourt les noeuds de la cellule actuelle
  841. C         °°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°
  842.  
  843. C         Pointeur de la liste des noeuds de la cellule actuelle
  844.           MLENTC = GRID3(IX0,IY0,IZ0)
  845. C         Nombre de noeuds contenus dans la cellule actuelle
  846.           NBPC0 = MLENTC.LECT(/1)
  847. C         Boucle sur chaque noeud de la cellule actuelle
  848.           DO 66, J = 1, NBPC0
  849. C             Récupération du noeud actuel
  850.               INDM0 = MLENTC.LECT(J)
  851.               NO0 = NUM(1,INDM0)
  852.               IND0 = (NO0-1)*(IDIM+1)
  853. C             Récupération des coordonnées du noeud actuel
  854.               X0 = XCOOR(IND0+1)
  855.               Y0 = XCOOR(IND0+2)
  856.               Z0 = XCOOR(IND0+3)
  857. C             Coordonnées -/+ R0
  858.               XLD0 = X0-R0
  859.               XLU0 = X0+R0
  860.               YLD0 = Y0-R0
  861.               YLU0 = Y0+R0
  862.               ZLD0 = Z0-R0
  863.               ZLU0 = Z0+R0
  864.  
  865. C             On génère la liste des noeuds voisins du noeud (MLENTN)
  866. C             du noeud actuel ainsi que leurs poids associés (MLREEP)
  867. C             -------------------------------------------------------
  868.  
  869. C             Initialisation de la liste des noeuds voisins (MLENTN)
  870. C             et de la liste de leurs poids respectifs (MLREEP)
  871.               JG = 1
  872.               SEGINI MLENTN, MLREEP
  873. 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)
  874.               MLENTN.LECT(1) = NO0
  875.               SELFW0 = VPOCHA(INDM0,1)
  876.               MLREEP.PROG(1) = SELFW0
  877. C             Initialisation de la somme des poids des noeuds voisins
  878.               SUMFAC0 = SELFW0
  879. C             On parcourt la liste des voisins potentiels
  880.               DO 34 K = 1, MLENTV.LECT(/1)
  881. C                 Récupération du noeud voisin
  882.                   INDM1 = MLENTV.LECT(K)
  883. C                 Si le noeud voisin traité est différent du noeud actuel
  884.                   IF(INDM1.NE.INDM0) THEN
  885.                       NO1 = NUM(1,INDM1)
  886.                       IND1 = (NO1-1)*(IDIM+1)
  887. C                     Coordonnées X du noeud voisin
  888.                       X1 = XCOOR(IND1+1)
  889. C                     Si la différence de coordonnées n'excède pas R0
  890.                       IF(X1.GE.XLD0.AND.X1.LE.XLU0) THEN
  891. C                         Coordonnées Y du noeud voisin
  892.                           Y1 = XCOOR(IND1+2)
  893. C                         Si la différence n'excède pas R0
  894.                           IF(Y1.GE.YLD0.AND.Y1.LE.YLU0) THEN
  895. C                             Coordonnées Z du noeud voisin
  896.                               Z1 = XCOOR(IND1+3)
  897. C                             Si la différence n'excède pas R0
  898.                               IF(Z1.GE.ZLD0.AND.Z1.LE.ZLU0) THEN
  899.                                   A0 = R0 - SQRT((X0-X1)**2 +
  900.      &                                       (Y0-Y1)**2 + (Z0-Z1)**2)
  901. C                                 Si la distance n'excède pas R0
  902.                                   IF(A0.GT.0.0) THEN
  903. C                                     Calcul du poids pour ce voisin
  904.                                       FAC0 = (A0/R0)**EMU0
  905.      &                                       * VPOCHA(INDM1,1)
  906. C                                     On sauvegarde le noeud voisin
  907. C                                     à la liste des noeuds voisins
  908. C                                     ainsi que son poids associé
  909.                                       JG = JG + 1
  910.                                       SEGADJ MLENTN, MLREEP
  911.                                       MLENTN.LECT(JG) = NO1
  912.                                       MLREEP.PROG(JG) = FAC0
  913. C                                     On met à jour la somme des poids
  914.                                       SUMFAC0 = SUMFAC0 + FAC0
  915.                                   ENDIF
  916.                               ENDIF
  917.                           ENDIF
  918.                       ENDIF
  919.                   ENDIF
  920. 34            CONTINUE
  921.  
  922. C             On ne conserve que les noeuds voisins de poids supérieur au critère
  923. C             -------------------------------------------------------------------
  924.  
  925. C             On boucle sur les noeuds voisins
  926.               A0 = CRI0 * SUMFAC0
  927.               DO 36 K = JG, 2, -1
  928.                   IF(MLREEP.PROG(K).LE.A0) THEN
  929. C                     On retire de la somme des poids
  930. C                     le poids du noeud voisin à éliminer
  931.                       SUMFAC0 = SUMFAC0 - MLREEP.PROG(K)
  932. C                     On élimine le noeud voisin et le poids associé
  933. C                      IF(K.NE.JG) THEN >> ce test est inutile car la boucle est ignorée si K>JG
  934.                           DO 38 L = K, JG-1
  935.                               MLENTN.LECT(L) = MLENTN.LECT(L+1)
  936.                               MLREEP.PROG(L) = MLREEP.PROG(L+1)
  937. 38                        CONTINUE
  938. C                      ENDIF
  939. C                     On met à jour le nombre de noeuds voisins
  940.                       JG = JG - 1
  941. C                     On ajuste la taille des listes
  942.                       SEGADJ MLENTN, MLREEP
  943.                   ENDIF
  944. 36            CONTINUE
  945.  
  946. C             Ajout d'une matrice élémentaire à la matrice de rigidité
  947. C             --------------------------------------------------------
  948.  
  949. C             On définit le nombre d'inconnues primales (= nombre de voisins)
  950.               NLIGRP = JG
  951. C             On définit le nombre de noeuds de l'élément de la matrice élémentaire
  952.               NBNN = JG
  953.  
  954.               IF(JG.GT.NBVMAX0) THEN
  955. C                 Si besoin on ajuste la taille de MLENTS
  956.                   NBVMAX0 = JG
  957.                   SEGADJ MLENTS
  958.               ENDIF
  959.  
  960. C             Les matrices élémentaires de même nombre de noeuds (ou voisins) sont regroupées dans une même sous-matrice.
  961. C             Si la sous-matrice pour ce nombre de voisins n'existe pas encore, on l'initialise
  962.               IF((JG.GT.NBVMAX0).OR.(MLENTS.LECT(JG).EQ.0)) THEN
  963. C                 On met à jour le nombre de sous matrices
  964.                   NRIGEL = NRIGEL + 1
  965. C                 Sauvegarde de l'indice de cette nouvelle sous-matrice
  966.                   MLENTS.LECT(JG) = NRIGEL
  967. C                 Initialisation du segment descripteur de la sous-matrice
  968. C                 On indique son nombre d'inconnues duales (déjà fait au début)
  969. C                  NLIGRD = 1
  970.                   SEGINI DESCR
  971. C                 La première et unique inconnue duale correspondant au 1er noeud voisin (le noeud actuel)
  972.                   NOELED(1) = 1
  973. C                 Nom de l'inconnue dual
  974.                   LISDUA(1) = NOMD0
  975. C                 Les inconnues primales, allant simplement de 1 au nombre de voisins,
  976. C                 portent toutes le même nom
  977.                   DO 39 L = 1, NLIGRP
  978.                       NOELEP(L) = L
  979.                       LISINC(L) = NOMP0
  980. 39                CONTINUE
  981. C                 Initialisation d'un maillage
  982.                   NBELEM = 1
  983.                   NBSOUS = 0
  984.                   NBREF = 0
  985.                   SEGINI IPT1
  986. C                 Initialisation de la sous-matrice
  987. C                 Elle contiendra au moins 1 élément, qui va être ajouté
  988.                   NELRIG = 1
  989.                    rigrel=0
  990.                   SEGINI XMATRI
  991. C                 On indique que la matrice ne possède pas de symétries
  992.                   SYMRE = 2
  993. C                 On ajuste le nombre de sous matrices
  994.                   JG2 = JG
  995.                   JG = NRIGEL
  996.                   SEGADJ MLENTD, MLENTM, MLENTX
  997.                   JG = JG2
  998. C                 Sauvegarde des pointeurs des 3 segments crées ci-dessus
  999.                   MLENTD.LECT(NRIGEL) = DESCR
  1000.                   MLENTM.LECT(NRIGEL) = IPT1
  1001.                   MLENTX.LECT(NRIGEL) = XMATRI
  1002. C             Si la sous-matrice existe, on la récupère
  1003.               ELSE
  1004.                   IPT1 = MLENTM.LECT(MLENTS.LECT(JG))
  1005.                   XMATRI = MLENTX.LECT(MLENTS.LECT(JG))
  1006. C                 On incrémente le nombre d'éléments du maillage
  1007. C                 et le nombre de matrices élémentaires de la sous-matrice
  1008.                   NBELEM = IPT1.NUM(/2) + 1
  1009.                   NELRIG = RE(/3) + 1
  1010.                    rigrel=0
  1011.                   SEGADJ IPT1, XMATRI
  1012.               ENDIF
  1013. C             On boucle sur les noeuds voisins
  1014.               DO 40 K = 1, JG
  1015. C                 Dans IPT1, on ajoute la liste de noeuds voisins au maillage
  1016.                   IPT1.NUM(K,NBELEM) = MLENTN.LECT(K)
  1017. C                 Dans XMATRI, on ajoute la liste des composantes de la nouvelle matrice de rigidité élémentaire
  1018.                   RE(1,K,NELRIG) = MLREEP.PROG(K)
  1019. 40            CONTINUE
  1020. 66        CONTINUE
  1021. 33    CONTINUE
  1022. C     On supprime les segments devenus inutiles
  1023.       SEGSUP GRID3SEG, INDEX3SEG
  1024.  
  1025. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1026. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1027. CC
  1028. CC  Combinaison des sous-matrices dans une matrice globale de rigidité
  1029. CC
  1030. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1031. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1032.  
  1033. 98    CONTINUE
  1034.  
  1035. C     On initialise le segment de la matrice de rigidité globale
  1036. C     qui prendra la taille de NRIGEL
  1037.       SEGINI MRIGID
  1038. C     On indique son type de matrice
  1039.       MTYMAT = 'RIGIDITE'
  1040.       IFORIG = mchpoi.IFOPOI
  1041.  
  1042. C     On remplit la matrice de rigidité avec les matrices élémentaires,
  1043. C     descripteurs et maillages correspondant à chaque sous-matrice
  1044.       DO 11 I = 1, NRIGEL
  1045. C               Récupération des segments du descripteur,
  1046. C               de la sous-matrice et de la liste des noeuds
  1047.                 DESCR  = MLENTD.LECT(I)
  1048.                 XMATRI = MLENTX.LECT(I)
  1049.                 IPT1   = MLENTM.LECT(I)
  1050. C               Affectation du maillage, du descripteur et du contenu de la sous-matrice
  1051.                 MRIGID.IRIGEL(1,I) = IPT1
  1052.                 MRIGID.IRIGEL(3,I) = DESCR
  1053.                 MRIGID.IRIGEL(4,I) = XMATRI
  1054. C               Retrait du caractère *MOD des SEGMENTS crees dans ce SP                
  1055.                 SEGACT,IPT1,DESCR,XMATRI
  1056. C               On indique que la matrice globale ne possède pas de symétries
  1057.                 MRIGID.IRIGEL(7,I) = 2
  1058. C               On précise son coefficient multiplicateur
  1059.                 MRIGID.COERIG(I) = 1.D0
  1060. 11    CONTINUE
  1061.       SEGACT,MRIGID
  1062.  
  1063. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1064. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1065. CC
  1066. CC  Nettoyage final et sortie
  1067. CC
  1068. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1069. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  1070.  
  1071. C     On supprime les segments utiles seulement au sein de la subroutine et quand R0 != 0
  1072.       SEGSUP MLENTX, MLENTD, MLENTM, MLENTN, MLREEP, NECX, NECY, NECZ
  1073.  
  1074. 99    CONTINUE
  1075.  
  1076.       SEGDES MCOORD
  1077.  
  1078. C     On désactive les segments restants avant de sortir
  1079.       SEGDES MRIGID, MELEME, MPOVAL, MSOUPO, MCHPOI
  1080. C     On renvoie la matrice de rigidité
  1081.       CALL ECROBJ('RIGIDITE',MRIGID)
  1082.  
  1083.       RETURN
  1084.       END
  1085.  
  1086.  
  1087.  

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