Télécharger la plaquette Cast3M


fr en

Télécharger dcov.eso

Retour à la liste

Numérotation des lignes : 

dcov

  1. C DCOV      SOURCE    CB215821  23/01/25    21:15:10     11573          
  2.       SUBROUTINE DCOV
  3. C
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8(A-H,O-Z)
  6. C
  7.  
  8. -INC PPARAM
  9. -INC CCOPTIO
  10. -INC SMELEME
  11. -INC SMCOORD
  12. -INC SMRIGID
  13. C
  14.       DIMENSION ZLAM(3),ZINVL(3),XCOV(3,3)
  15.       CHARACTER*8 MOTS(1),MOTLA(4)
  16.       CHARACTER*4 MOTLO(2)
  17.       CHARACTER*9 MOTD(1)
  18. C
  19.       DATA MOTS/'SIGMA   '/
  20.       DATA MOTLA/'LAMBDA  ','LAMBDA1 ','LAMBDA2 ','LAMBDA3 '/
  21.       DATA MOTLO/'EXPO','GAUS'/
  22.       DATA MOTD/'DIRECTION'/
  23. C
  24. C  epsilon servant à différents tests
  25. C
  26.       EPS = 1.D-12
  27. C
  28.        IDIR=0
  29. C
  30. C  Lecture du maillage
  31. C
  32.       CALL LIROBJ('MAILLAGE',MELEME,1,IRET)
  33.       IF (IERR.NE.0) RETURN
  34. C
  35. C  Lecture de la loi utilisée
  36. C
  37.       CALL LIRMOT(MOTLO,2,ILOI,1)
  38.       IF (ILOI.EQ.0) RETURN
  39. C
  40. C  Lecture du mot-clé 'SIGMA'
  41. C
  42.       CALL LIRMOT(MOTS,1,IRET,1)
  43.       IF (IRET.EQ.0) RETURN
  44. C
  45. C  Lecture de la valeur de sigma (strictement supérieure à 0.D0)
  46. C
  47.       CALL LIRREE(ZSIG,1,IRET)
  48.       IF (IERR.NE.0) RETURN
  49.       IF (ZSIG.LE.0.D0) THEN
  50.         REAERR(1) = REAL(ZSIG)
  51.         REAERR(2) = REAL(0.D0)
  52.         MOTERR(1:8) = ' SIGMA  '
  53.         CALL ERREUR(41)
  54.         RETURN
  55.       ENDIF
  56. C
  57. C  Lecture du mot-clé 'LAMBDA' ou 'LAMBDA1'
  58. C
  59.       IF (IDIM.EQ.1) THEN
  60.         CALL LIRMOT(MOTLA,1,IANI,1)
  61.       ELSE
  62.         CALL LIRMOT(MOTLA,2,IANI,1)
  63.       ENDIF
  64.       IF (IANI.EQ.0) RETURN
  65. C
  66. C  Lecture de la valeur de lambda (ou lambda1)
  67. C                                        (strictement supérieure à 0.D0)
  68. C
  69.       CALL LIRREE(ZLAM(1),1,IRET)
  70.       IF (IERR.NE.0) RETURN
  71.       IF (ZLAM(1).LE.0.D0) THEN
  72.         REAERR(1) = REAL(ZLAM(1))
  73.         REAERR(2) = REAL(0.D0)
  74.         IF (IANI.EQ.1) MOTERR(1:8) = ' LAMBDA '
  75.         IF (IANI.EQ.2) MOTERR(1:8) = ' LAMBDA1'
  76.         CALL ERREUR(41)
  77.         RETURN
  78.       ENDIF
  79. C
  80.  
  81.       IF (IANI.EQ.2) THEN
  82. C
  83. C  Lecture éventuelle du mot-clé 'LAMBDA2' (cas anisotrope)
  84. C
  85.         CALL LIRMOT(MOTLA(3),1,IRET,1)
  86.         IF (IRET.EQ.0) RETURN
  87. C
  88. C  Lecture éventuelle de la valeur de lambda2 (cas anisotrope)
  89. C                                        (strictement supérieure à 0.D0)
  90. C
  91.         CALL LIRREE(ZLAM(2),1,IRET)
  92.         IF (IERR.NE.0) RETURN
  93.         IF (ZLAM(2).LE.0.D0) THEN
  94.           REAERR(1) = REAL(ZLAM(2))
  95.           REAERR(2) = REAL(0.D0)
  96.           MOTERR(1:8) = ' LAMBDA2'
  97.           CALL ERREUR(41)
  98.           RETURN
  99.         ENDIF
  100. C
  101.         IF (IDIM.EQ.3) THEN
  102. C
  103. C  Lecture éventuelle du mot-clé 'LAMBDA3' (cas anisotrope en 3d)
  104. C
  105.           CALL LIRMOT(MOTLA(4),1,IRET,1)
  106.           IF (IRET.EQ.0) RETURN
  107. C
  108. C  Lecture éventuelle de la valeur de lambda3 (cas anisotrope en 3d)
  109. C                                        (strictement supérieure à 0.D0)
  110. C
  111.           CALL LIRREE(ZLAM(3),1,IRET)
  112.           IF (IERR.NE.0) RETURN
  113.           IF (ZLAM(3).LE.0.D0) THEN
  114.             REAERR(1) = REAL(ZLAM(3))
  115.             REAERR(2) = REAL(0.D0)
  116.             MOTERR(1:8) = ' LAMBDA3'
  117.             CALL ERREUR(41)
  118.             RETURN
  119.           ENDIF
  120. C
  121.         ENDIF
  122. C
  123. C  Lecture éventuelle et optionnelle du mot-clé 'DIRECTION'
  124. C     (cas anisotrope, directions d'anisotropie non parallèles aux axes)
  125. C
  126.         CALL LIRMOT(MOTD,1,IDIR,0)
  127.         IF (IERR.NE.0) RETURN
  128. C
  129.         IF (IDIR.NE.0) THEN
  130. C
  131. C  Lecture éventuelle de la valeur de VEC1
  132. C     (cas anisotrope, directions d'anisotropie non parallèles aux axes)
  133. C
  134.           CALL LIROBJ('POINT',IPTV1,1,IRET)
  135.           IF (IERR.NE.0) RETURN
  136. C
  137.           IF (IDIM.EQ.3) THEN
  138. C
  139. C  Lecture éventuelle de la valeur de VEC2
  140. C     (cas anisotrope en 3d,
  141. C                      directions d'anisotropie non parallèles aux axes)
  142. C
  143.             CALL LIROBJ('POINT',IPTV2,1,IRET)
  144.             IF (IERR.NE.0) RETURN
  145. C
  146.           ENDIF
  147. C
  148.         ENDIF
  149. C
  150.       ENDIF
  151. C
  152. C  Dans le cas anisotrope,
  153. C                      directions d'anisotropie non parallèles aux axes,
  154. C  on vérifie que le(s) vecteur(s) introduit(s) sous forme de point ne
  155. C  sont pas de longueur nulle
  156. C  si tel est le cas, on le(s) norme
  157. C
  158.       SEGACT,MCOORD
  159.       IF (IDIR.NE.0) THEN
  160. C
  161.         IREF1 = (IPTV1-1)*(IDIM+1)
  162.         SCOV1 = 0.D0
  163.         DO 1 I=1,IDIM
  164.           XCOV(1,I) = XCOOR(IREF1+I)
  165.           SCOV1 = SCOV1+(XCOV(1,I)*XCOV(1,I))
  166.     1   CONTINUE
  167. C
  168.         IF (SCOV1.LT.EPS) THEN
  169.           CALL ERREUR(239)
  170.           RETURN
  171.         ELSE
  172.           DO 2 I=1,IDIM
  173.             XCOV(1,I) = XCOV(1,I)/SCOV1
  174.     2     CONTINUE
  175.         ENDIF
  176. C
  177.         IF (IDIM.EQ.3) THEN
  178. C
  179.           IREF2 = (IPTV2-1)*(IDIM+1)
  180.           SCOV2 = 0.D0
  181.           DO 3 I=1,IDIM
  182.             XCOV(2,I) = XCOOR(IREF2+I)
  183.             SCOV2 = SCOV2+(XCOV(2,I)*XCOV(2,I))
  184.     3     CONTINUE
  185. C
  186.           IF (SCOV2.LT.EPS) THEN
  187.             CALL ERREUR(239)
  188.             RETURN
  189.           ELSE
  190.             DO 4 I=1,IDIM
  191.               XCOV(2,I) = XCOV(2,I)/SCOV2
  192.     4       CONTINUE
  193.           ENDIF
  194. C
  195. C  Dans le cas anisotrope 3d,
  196. C                      directions d'anisotropie non parallèles aux axes,
  197. C  on vérifie que les vecteurs introduits sous forme de point
  198. C  sont bien orthogonaux
  199. C
  200.           SCOV12 = 0.D0
  201.           DO 5 I=1,IDIM
  202.             SCOV12 = SCOV12+(XCOV(1,I)*XCOV(2,I))
  203.     5     CONTINUE
  204.           IF (SCOV12.GT.EPS) THEN
  205.             CALL ERREUR(161)
  206.             RETURN
  207.           ENDIF
  208. C
  209.         ENDIF
  210. C
  211.       ENDIF
  212. C
  213.       IF (IDIR.NE.0) THEN
  214. C
  215. C  Dans le cas anisotrope 2d,
  216. C                      directions d'anisotropie non parallèles aux axes,
  217. C  on complète le repère par rotation de +90 degrés
  218. C
  219.         IF (IDIM.EQ.2) THEN
  220.           XCOV(2,1) = -1.D0*XCOV(1,2)
  221.           XCOV(2,2) = XCOV(1,1)
  222.         ENDIF
  223. C
  224. C  Dans le cas anisotrope 3d,
  225. C                      directions d'anisotropie non parallèles aux axes,
  226. C  on complète le repère par produit vectoriel des deux premiers
  227. C  vecteurs
  228. C
  229.         IF (IDIM.EQ.3) THEN
  230.           XCOV(3,1) = (XCOV(1,2)*XCOV(2,3))-(XCOV(1,3)*XCOV(2,2))
  231.           XCOV(3,2) = (XCOV(1,3)*XCOV(2,1))-(XCOV(1,1)*XCOV(2,3))
  232.           XCOV(3,3) = (XCOV(1,1)*XCOV(2,2))-(XCOV(1,2)*XCOV(2,1))
  233.         ENDIF
  234. C
  235.       ENDIF
  236. C
  237. C  Il vaut mieux faire * (1./lambda) que / lambda
  238. C
  239.       ZINVL(1) = 1.D0/ZLAM(1)
  240.       IF (IANI.EQ.2) THEN
  241.         ZINVL(2) = 1.D0/ZLAM(2)
  242.         IF (IDIM.EQ.3) ZINVL(3) = 1.D0/ZLAM(3)
  243.       ENDIF
  244. C
  245. C  Les éléments du MELEME sont transformés en POI1 si nécessaire
  246. C
  247.       SEGACT MELEME
  248.       IPT1=MELEME
  249.       IPT2=MELEME
  250.       IF(ITYPEL.NE.1.OR.(LISOUS(/1)).NE.0) THEN
  251.         CALL CHANGE(IPT1,1)
  252.         IF (IERR.NE.0) RETURN
  253.       ENDIF
  254. C
  255.       NBLI=IPT1.NUM(/2)
  256. C
  257. C  On définit un super-élément correspondant aux noeuds
  258. C
  259.       NBSOUS = 0
  260.       NBREF = 0
  261.       NBNN = NBLI
  262.       NBELEM = 1
  263.       SEGINI MELEME
  264.       ITYPEL = 28
  265.       DO 7 I=1,NBLI
  266.         NUM(I,1) = IPT1.NUM(1,I)
  267.     7 CONTINUE
  268.       SEGDES IPT2
  269. C
  270. C  On met la matrice dans un objet rigidité
  271. C
  272.       NRIGE = 8
  273.       NRIGEL = 1
  274.       SEGINI MRIGID
  275.       ICHOLE = 0
  276.       IMGEO1 = 0
  277.       IMGEO2 = 0
  278.       IFORIG = IFOUR
  279.       ISUPEQ = 0
  280.       NELRIG = 1
  281.       COERIG(1) = 1.D0
  282. *      SEGINI IMATRI
  283.       NLIGRD = NBLI
  284.       NLIGRP = NBLI
  285.           segini xmatri
  286.       SEGINI DESCR
  287.       IRIGEL(1,1) = MELEME
  288.       IRIGEL(2,1) = 0
  289.       IRIGEL(3,1) = DESCR
  290.       IRIGEL(4,1) = xMATRI
  291.       IRIGEL(5,1) = 0
  292.       IRIGEL(6,1) = 0
  293.       IRIGEL(7,1) = 2
  294.       DO 8 I=1,NBLI
  295.         NOELEP(I) = I
  296.         NOELED(I) = I
  297.         LISINC(I) = 'SCAL'
  298.         LISDUA(I) = 'SCAL'
  299.     8 CONTINUE
  300.       SEGDES DESCR,MRIGID
  301. *      SEGINI XMATRI
  302. *      IMATTT(1) = XMATRI
  303. *      DO 9 I=1,NBLI
  304. *        DO 10 J=1,NBLI
  305. *          RE(I,J,1) = 0.D0
  306. *   10   CONTINUE
  307. *    9 CONTINUE
  308. C
  309.       IF (ILOI.EQ.1) THEN
  310. C
  311. C Loi exponentielle
  312. C
  313.         IF (IANI.EQ.1) THEN
  314. C
  315. C Cas isotrope
  316. C
  317.           DO 100 I=1,NBLI
  318.             IREF=(NUM(I,1)-1)*(IDIM+1)
  319.             DO 110 J=1,I
  320.               JREF=(NUM(J,1)-1)*(IDIM+1)
  321.               VM=0.D0
  322.               DO 120 I1=1,IDIM
  323.                 VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2
  324.   120         CONTINUE
  325.               VM=SQRT(VM)
  326.               RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1))
  327.   110       CONTINUE
  328.   100     CONTINUE
  329. C
  330.         ELSEIF (IANI.EQ.2) THEN
  331. C
  332. C Cas anisotrope :
  333. C  directions d'anisotropie parallèles aux axes (IDIR nul)
  334. C
  335.           IF (IDIR.EQ.0) THEN
  336. C
  337.             DO 200 I=1,NBLI
  338.               IREF=(NUM(I,1)-1)*(IDIM+1)
  339.               DO 210 J=1,I
  340.                 JREF=(NUM(J,1)-1)*(IDIM+1)
  341.                 VM=0.D0
  342.                 DO 220 I1=1,IDIM
  343.                   VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2
  344.   220           CONTINUE
  345.                 VM=SQRT(VM)
  346.                 RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  347.   210         CONTINUE
  348.   200       CONTINUE
  349. C
  350. C Cas anisotrope :
  351. C  directions d'anisotropie non parallèles aux axes (IDIR non nul)
  352. C
  353.           ELSEIF (IDIR.NE.0) THEN
  354. C
  355.             DO 300 I=1,NBLI
  356.               IREF=(NUM(I,1)-1)*(IDIM+1)
  357.               DO 310 J=1,I
  358.                 JREF=(NUM(J,1)-1)*(IDIM+1)
  359.                 VM=0.D0
  360.                 DO 320 I1=1,IDIM
  361.                   PSCAL=0.D0
  362.                   DO 330 J1=1,IDIM
  363.                     PSCAL = PSCAL +
  364.      *                       (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1)
  365.   330             CONTINUE
  366.                   VM=VM+(PSCAL*ZINVL(I1))**2
  367.   320           CONTINUE
  368.                 VM=SQRT(VM)
  369.                 RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  370.   310         CONTINUE
  371.   300       CONTINUE
  372. C
  373.           ENDIF
  374. C
  375.         ENDIF
  376. C
  377.       ELSEIF (ILOI.EQ.2) THEN
  378. C
  379. C Loi gaussienne
  380. C
  381.         IF (IANI.EQ.1) THEN
  382. C
  383. C Cas isotrope
  384. C
  385.           DO 400 I=1,NBLI
  386.             IREF=(NUM(I,1)-1)*(IDIM+1)
  387.             DO 410 J=1,I
  388.               JREF=(NUM(J,1)-1)*(IDIM+1)
  389.               VM=0.D0
  390.               DO 420 I1=1,IDIM
  391.                 VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2
  392.   420         CONTINUE
  393.               RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1))
  394.   410       CONTINUE
  395.   400     CONTINUE
  396. C
  397.         ELSEIF (IANI.EQ.2) THEN
  398. C
  399. C Cas anisotrope :
  400. C  directions d'anisotropie parallèles aux axes (IDIR nul)
  401. C
  402.           IF (IDIR.EQ.0) THEN
  403. C
  404.             DO 500 I=1,NBLI
  405.               IREF=(NUM(I,1)-1)*(IDIM+1)
  406.               DO 510 J=1,I
  407.                 JREF=(NUM(J,1)-1)*(IDIM+1)
  408.                 VM=0.D0
  409.                 DO 520 I1=1,IDIM
  410.                   VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2
  411.   520           CONTINUE
  412.                 RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  413.   510         CONTINUE
  414.   500       CONTINUE
  415. C
  416. C Cas anisotrope :
  417. C  directions d'anisotropie non parallèles aux axes (IDIR non nul)
  418. C
  419.           ELSEIF (IDIR.NE.0) THEN
  420. C
  421.             DO 600 I=1,NBLI
  422.               IREF=(NUM(I,1)-1)*(IDIM+1)
  423.               DO 610 J=1,I
  424.                 JREF=(NUM(J,1)-1)*(IDIM+1)
  425.                 VM=0.D0
  426.                 DO 620 I1=1,IDIM
  427.                   PSCAL=0.D0
  428.                   DO 630 J1=1,IDIM
  429.                     PSCAL = PSCAL +
  430.      *                       (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1)
  431.   630             CONTINUE
  432.                   VM=VM+(PSCAL*ZINVL(I1))**2
  433.   620           CONTINUE
  434.                 RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  435.   610         CONTINUE
  436.   600       CONTINUE
  437. C
  438.           ENDIF
  439. C
  440.         ENDIF
  441. C
  442.       ENDIF
  443.       SEGDES,MCOORD
  444. C
  445. C Calcul de M telle que (RE =) COV = M Mt
  446. C   M est stockée dans la partie inférieure de RE
  447. C
  448.       CALL SYCHOM(RE,NBLI)
  449. C
  450. C Création de l'objet rigidité
  451. C
  452.       CALL ACTOBJ('RIGIDITE',MRIGID,0)
  453.       CALL ECROBJ('RIGIDITE',MRIGID)
  454. C
  455.       END
  456.  
  457.  

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