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

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