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    PV090527  26/04/30    21:15:27     12529          
  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.           rigrel=0
  286.           segini xmatri
  287.       SEGINI DESCR
  288.       IRIGEL(1,1) = MELEME
  289.       IRIGEL(2,1) = 0
  290.       IRIGEL(3,1) = DESCR
  291.       IRIGEL(4,1) = xMATRI
  292.       IRIGEL(5,1) = 0
  293.       IRIGEL(6,1) = 0
  294.       IRIGEL(7,1) = 2
  295.       DO 8 I=1,NBLI
  296.         NOELEP(I) = I
  297.         NOELED(I) = I
  298.         LISINC(I) = 'SCAL'
  299.         LISDUA(I) = 'SCAL'
  300.     8 CONTINUE
  301.       SEGDES DESCR,MRIGID
  302. *      SEGINI XMATRI
  303. *      IMATTT(1) = XMATRI
  304. *      DO 9 I=1,NBLI
  305. *        DO 10 J=1,NBLI
  306. *          RE(I,J,1) = 0.D0
  307. *   10   CONTINUE
  308. *    9 CONTINUE
  309. C
  310.       IF (ILOI.EQ.1) THEN
  311. C
  312. C Loi exponentielle
  313. C
  314.         IF (IANI.EQ.1) THEN
  315. C
  316. C Cas isotrope
  317. C
  318.           DO 100 I=1,NBLI
  319.             IREF=(NUM(I,1)-1)*(IDIM+1)
  320.             DO 110 J=1,I
  321.               JREF=(NUM(J,1)-1)*(IDIM+1)
  322.               VM=0.D0
  323.               DO 120 I1=1,IDIM
  324.                 VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2
  325.   120         CONTINUE
  326.               VM=SQRT(VM)
  327.               RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1))
  328.   110       CONTINUE
  329.   100     CONTINUE
  330. C
  331.         ELSEIF (IANI.EQ.2) THEN
  332. C
  333. C Cas anisotrope :
  334. C  directions d'anisotropie parallèles aux axes (IDIR nul)
  335. C
  336.           IF (IDIR.EQ.0) THEN
  337. C
  338.             DO 200 I=1,NBLI
  339.               IREF=(NUM(I,1)-1)*(IDIM+1)
  340.               DO 210 J=1,I
  341.                 JREF=(NUM(J,1)-1)*(IDIM+1)
  342.                 VM=0.D0
  343.                 DO 220 I1=1,IDIM
  344.                   VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2
  345.   220           CONTINUE
  346.                 VM=SQRT(VM)
  347.                 RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  348.   210         CONTINUE
  349.   200       CONTINUE
  350. C
  351. C Cas anisotrope :
  352. C  directions d'anisotropie non parallèles aux axes (IDIR non nul)
  353. C
  354.           ELSEIF (IDIR.NE.0) THEN
  355. C
  356.             DO 300 I=1,NBLI
  357.               IREF=(NUM(I,1)-1)*(IDIM+1)
  358.               DO 310 J=1,I
  359.                 JREF=(NUM(J,1)-1)*(IDIM+1)
  360.                 VM=0.D0
  361.                 DO 320 I1=1,IDIM
  362.                   PSCAL=0.D0
  363.                   DO 330 J1=1,IDIM
  364.                     PSCAL = PSCAL +
  365.      *                       (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1)
  366.   330             CONTINUE
  367.                   VM=VM+(PSCAL*ZINVL(I1))**2
  368.   320           CONTINUE
  369.                 VM=SQRT(VM)
  370.                 RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  371.   310         CONTINUE
  372.   300       CONTINUE
  373. C
  374.           ENDIF
  375. C
  376.         ENDIF
  377. C
  378.       ELSEIF (ILOI.EQ.2) THEN
  379. C
  380. C Loi gaussienne
  381. C
  382.         IF (IANI.EQ.1) THEN
  383. C
  384. C Cas isotrope
  385. C
  386.           DO 400 I=1,NBLI
  387.             IREF=(NUM(I,1)-1)*(IDIM+1)
  388.             DO 410 J=1,I
  389.               JREF=(NUM(J,1)-1)*(IDIM+1)
  390.               VM=0.D0
  391.               DO 420 I1=1,IDIM
  392.                 VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2
  393.   420         CONTINUE
  394.               RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1))
  395.   410       CONTINUE
  396.   400     CONTINUE
  397. C
  398.         ELSEIF (IANI.EQ.2) THEN
  399. C
  400. C Cas anisotrope :
  401. C  directions d'anisotropie parallèles aux axes (IDIR nul)
  402. C
  403.           IF (IDIR.EQ.0) THEN
  404. C
  405.             DO 500 I=1,NBLI
  406.               IREF=(NUM(I,1)-1)*(IDIM+1)
  407.               DO 510 J=1,I
  408.                 JREF=(NUM(J,1)-1)*(IDIM+1)
  409.                 VM=0.D0
  410.                 DO 520 I1=1,IDIM
  411.                   VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2
  412.   520           CONTINUE
  413.                 RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  414.   510         CONTINUE
  415.   500       CONTINUE
  416. C
  417. C Cas anisotrope :
  418. C  directions d'anisotropie non parallèles aux axes (IDIR non nul)
  419. C
  420.           ELSEIF (IDIR.NE.0) THEN
  421. C
  422.             DO 600 I=1,NBLI
  423.               IREF=(NUM(I,1)-1)*(IDIM+1)
  424.               DO 610 J=1,I
  425.                 JREF=(NUM(J,1)-1)*(IDIM+1)
  426.                 VM=0.D0
  427.                 DO 620 I1=1,IDIM
  428.                   PSCAL=0.D0
  429.                   DO 630 J1=1,IDIM
  430.                     PSCAL = PSCAL +
  431.      *                       (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1)
  432.   630             CONTINUE
  433.                   VM=VM+(PSCAL*ZINVL(I1))**2
  434.   620           CONTINUE
  435.                 RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  436.   610         CONTINUE
  437.   600       CONTINUE
  438. C
  439.           ENDIF
  440. C
  441.         ENDIF
  442. C
  443.       ENDIF
  444.       SEGDES,MCOORD
  445. C
  446. C Calcul de M telle que (RE =) COV = M Mt
  447. C   M est stockée dans la partie inférieure de RE
  448. C
  449.       CALL SYCHOM(RE,NBLI)
  450. C
  451. C Création de l'objet rigidité
  452. C
  453.       CALL ACTOBJ('RIGIDITE',MRIGID,0)
  454.       CALL ECROBJ('RIGIDITE',MRIGID)
  455. C
  456.       END
  457.  
  458.  
  459.  
  460.  

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