C DCOV      SOURCE    PV090527  26/04/30    21:15:27     12529          
      SUBROUTINE DCOV
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C

-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMCOORD
-INC SMRIGID
C
      DIMENSION ZLAM(3),ZINVL(3),XCOV(3,3)
      CHARACTER*8 MOTS(1),MOTLA(4)
      CHARACTER*4 MOTLO(2)
      CHARACTER*9 MOTD(1)
C
      DATA MOTS/'SIGMA   '/
      DATA MOTLA/'LAMBDA  ','LAMBDA1 ','LAMBDA2 ','LAMBDA3 '/
      DATA MOTLO/'EXPO','GAUS'/
      DATA MOTD/'DIRECTION'/
C
C  epsilon servant à différents tests
C
      EPS = 1.D-12
C
       IDIR=0
C
C  Lecture du maillage
C
      CALL LIROBJ('MAILLAGE',MELEME,1,IRET)
      IF (IERR.NE.0) RETURN
C
C  Lecture de la loi utilisée
C
      CALL LIRMOT(MOTLO,2,ILOI,1)
      IF (ILOI.EQ.0) RETURN
C
C  Lecture du mot-clé 'SIGMA'
C
      CALL LIRMOT(MOTS,1,IRET,1)
      IF (IRET.EQ.0) RETURN
C
C  Lecture de la valeur de sigma (strictement supérieure à 0.D0)
C
      CALL LIRREE(ZSIG,1,IRET)
      IF (IERR.NE.0) RETURN
      IF (ZSIG.LE.0.D0) THEN
        REAERR(1) = REAL(ZSIG)
        REAERR(2) = REAL(0.D0)
        MOTERR(1:8) = ' SIGMA  '
        CALL ERREUR(41)
        RETURN
      ENDIF
C
C  Lecture du mot-clé 'LAMBDA' ou 'LAMBDA1'
C
      IF (IDIM.EQ.1) THEN
        CALL LIRMOT(MOTLA,1,IANI,1)
      ELSE
        CALL LIRMOT(MOTLA,2,IANI,1)
      ENDIF
      IF (IANI.EQ.0) RETURN
C
C  Lecture de la valeur de lambda (ou lambda1)
C                                        (strictement supérieure à 0.D0)
C
      CALL LIRREE(ZLAM(1),1,IRET)
      IF (IERR.NE.0) RETURN
      IF (ZLAM(1).LE.0.D0) THEN
        REAERR(1) = REAL(ZLAM(1))
        REAERR(2) = REAL(0.D0)
        IF (IANI.EQ.1) MOTERR(1:8) = ' LAMBDA '
        IF (IANI.EQ.2) MOTERR(1:8) = ' LAMBDA1'
        CALL ERREUR(41)
        RETURN
      ENDIF
C
      
      IF (IANI.EQ.2) THEN
C
C  Lecture éventuelle du mot-clé 'LAMBDA2' (cas anisotrope)
C
        CALL LIRMOT(MOTLA(3),1,IRET,1)
        IF (IRET.EQ.0) RETURN
C
C  Lecture éventuelle de la valeur de lambda2 (cas anisotrope)
C                                        (strictement supérieure à 0.D0)
C
        CALL LIRREE(ZLAM(2),1,IRET)
        IF (IERR.NE.0) RETURN
        IF (ZLAM(2).LE.0.D0) THEN
          REAERR(1) = REAL(ZLAM(2))
          REAERR(2) = REAL(0.D0)
          MOTERR(1:8) = ' LAMBDA2'
          CALL ERREUR(41)
          RETURN
        ENDIF
C
        IF (IDIM.EQ.3) THEN
C
C  Lecture éventuelle du mot-clé 'LAMBDA3' (cas anisotrope en 3d)
C
          CALL LIRMOT(MOTLA(4),1,IRET,1)
          IF (IRET.EQ.0) RETURN
C
C  Lecture éventuelle de la valeur de lambda3 (cas anisotrope en 3d)
C                                        (strictement supérieure à 0.D0)
C
          CALL LIRREE(ZLAM(3),1,IRET)
          IF (IERR.NE.0) RETURN
          IF (ZLAM(3).LE.0.D0) THEN
            REAERR(1) = REAL(ZLAM(3))
            REAERR(2) = REAL(0.D0)
            MOTERR(1:8) = ' LAMBDA3'
            CALL ERREUR(41)
            RETURN
          ENDIF
C
        ENDIF
C
C  Lecture éventuelle et optionnelle du mot-clé 'DIRECTION'
C     (cas anisotrope, directions d'anisotropie non parallèles aux axes)
C
        CALL LIRMOT(MOTD,1,IDIR,0)
        IF (IERR.NE.0) RETURN
C
        IF (IDIR.NE.0) THEN
C
C  Lecture éventuelle de la valeur de VEC1
C     (cas anisotrope, directions d'anisotropie non parallèles aux axes)
C
          CALL LIROBJ('POINT',IPTV1,1,IRET)
          IF (IERR.NE.0) RETURN
C
          IF (IDIM.EQ.3) THEN
C
C  Lecture éventuelle de la valeur de VEC2
C     (cas anisotrope en 3d,
C                      directions d'anisotropie non parallèles aux axes)
C
            CALL LIROBJ('POINT',IPTV2,1,IRET)
            IF (IERR.NE.0) RETURN
C
          ENDIF
C
        ENDIF
C
      ENDIF
C
C  Dans le cas anisotrope,
C                      directions d'anisotropie non parallèles aux axes,
C  on vérifie que le(s) vecteur(s) introduit(s) sous forme de point ne
C  sont pas de longueur nulle
C  si tel est le cas, on le(s) norme
C
      SEGACT,MCOORD
      IF (IDIR.NE.0) THEN
C
        IREF1 = (IPTV1-1)*(IDIM+1)
        SCOV1 = 0.D0
        DO 1 I=1,IDIM
          XCOV(1,I) = XCOOR(IREF1+I)
          SCOV1 = SCOV1+(XCOV(1,I)*XCOV(1,I))
    1   CONTINUE
C
        IF (SCOV1.LT.EPS) THEN
          CALL ERREUR(239)
          RETURN
        ELSE
          DO 2 I=1,IDIM
            XCOV(1,I) = XCOV(1,I)/SCOV1
    2     CONTINUE
        ENDIF
C
        IF (IDIM.EQ.3) THEN
C
          IREF2 = (IPTV2-1)*(IDIM+1)
          SCOV2 = 0.D0
          DO 3 I=1,IDIM
            XCOV(2,I) = XCOOR(IREF2+I)
            SCOV2 = SCOV2+(XCOV(2,I)*XCOV(2,I))
    3     CONTINUE
C
          IF (SCOV2.LT.EPS) THEN
            CALL ERREUR(239)
            RETURN
          ELSE
            DO 4 I=1,IDIM
              XCOV(2,I) = XCOV(2,I)/SCOV2
    4       CONTINUE
          ENDIF
C
C  Dans le cas anisotrope 3d,
C                      directions d'anisotropie non parallèles aux axes,
C  on vérifie que les vecteurs introduits sous forme de point
C  sont bien orthogonaux
C
          SCOV12 = 0.D0
          DO 5 I=1,IDIM
            SCOV12 = SCOV12+(XCOV(1,I)*XCOV(2,I))
    5     CONTINUE
          IF (SCOV12.GT.EPS) THEN
            CALL ERREUR(161)
            RETURN
          ENDIF
C
        ENDIF
C
      ENDIF
C
      IF (IDIR.NE.0) THEN
C
C  Dans le cas anisotrope 2d,
C                      directions d'anisotropie non parallèles aux axes,
C  on complète le repère par rotation de +90 degrés
C
        IF (IDIM.EQ.2) THEN
          XCOV(2,1) = -1.D0*XCOV(1,2)
          XCOV(2,2) = XCOV(1,1)
        ENDIF
C
C  Dans le cas anisotrope 3d,
C                      directions d'anisotropie non parallèles aux axes,
C  on complète le repère par produit vectoriel des deux premiers
C  vecteurs
C
        IF (IDIM.EQ.3) THEN
          XCOV(3,1) = (XCOV(1,2)*XCOV(2,3))-(XCOV(1,3)*XCOV(2,2))
          XCOV(3,2) = (XCOV(1,3)*XCOV(2,1))-(XCOV(1,1)*XCOV(2,3))
          XCOV(3,3) = (XCOV(1,1)*XCOV(2,2))-(XCOV(1,2)*XCOV(2,1))
        ENDIF
C
      ENDIF
C
C  Il vaut mieux faire * (1./lambda) que / lambda
C
      ZINVL(1) = 1.D0/ZLAM(1)
      IF (IANI.EQ.2) THEN
        ZINVL(2) = 1.D0/ZLAM(2)
        IF (IDIM.EQ.3) ZINVL(3) = 1.D0/ZLAM(3)
      ENDIF
C
C  Les éléments du MELEME sont transformés en POI1 si nécessaire
C
      SEGACT MELEME
      IPT1=MELEME
      IPT2=MELEME
      IF(ITYPEL.NE.1.OR.(LISOUS(/1)).NE.0) THEN
        CALL CHANGE(IPT1,1)
        IF (IERR.NE.0) RETURN
      ENDIF
C
      NBLI=IPT1.NUM(/2)
C
C  On définit un super-élément correspondant aux noeuds
C
      NBSOUS = 0
      NBREF = 0
      NBNN = NBLI
      NBELEM = 1
      SEGINI MELEME
      ITYPEL = 28
      DO 7 I=1,NBLI
        NUM(I,1) = IPT1.NUM(1,I)
    7 CONTINUE
      SEGDES IPT2
C
C  On met la matrice dans un objet rigidité
C
      NRIGE = 8
      NRIGEL = 1
      SEGINI MRIGID
      ICHOLE = 0
      IMGEO1 = 0
      IMGEO2 = 0
      IFORIG = IFOUR
      ISUPEQ = 0
      NELRIG = 1
      COERIG(1) = 1.D0
*      SEGINI IMATRI
      NLIGRD = NBLI
      NLIGRP = NBLI
          rigrel=0
          segini xmatri
      SEGINI DESCR
      IRIGEL(1,1) = MELEME
      IRIGEL(2,1) = 0
      IRIGEL(3,1) = DESCR
      IRIGEL(4,1) = xMATRI
      IRIGEL(5,1) = 0
      IRIGEL(6,1) = 0
      IRIGEL(7,1) = 2
      DO 8 I=1,NBLI
        NOELEP(I) = I
        NOELED(I) = I
        LISINC(I) = 'SCAL'
        LISDUA(I) = 'SCAL'
    8 CONTINUE
      SEGDES DESCR,MRIGID
*      SEGINI XMATRI
*      IMATTT(1) = XMATRI
*      DO 9 I=1,NBLI
*        DO 10 J=1,NBLI
*          RE(I,J,1) = 0.D0
*   10   CONTINUE
*    9 CONTINUE
C
      IF (ILOI.EQ.1) THEN
C
C Loi exponentielle
C
        IF (IANI.EQ.1) THEN
C
C Cas isotrope
C
          DO 100 I=1,NBLI
            IREF=(NUM(I,1)-1)*(IDIM+1)
            DO 110 J=1,I
              JREF=(NUM(J,1)-1)*(IDIM+1)
              VM=0.D0
              DO 120 I1=1,IDIM
                VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2
  120         CONTINUE
              VM=SQRT(VM)
              RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1))
  110       CONTINUE
  100     CONTINUE
C
        ELSEIF (IANI.EQ.2) THEN
C
C Cas anisotrope :
C  directions d'anisotropie parallèles aux axes (IDIR nul)
C
          IF (IDIR.EQ.0) THEN
C
            DO 200 I=1,NBLI
              IREF=(NUM(I,1)-1)*(IDIM+1)
              DO 210 J=1,I
                JREF=(NUM(J,1)-1)*(IDIM+1)
                VM=0.D0
                DO 220 I1=1,IDIM
                  VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2
  220           CONTINUE
                VM=SQRT(VM)
                RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  210         CONTINUE
  200       CONTINUE
C
C Cas anisotrope :
C  directions d'anisotropie non parallèles aux axes (IDIR non nul)
C
          ELSEIF (IDIR.NE.0) THEN
C
            DO 300 I=1,NBLI
              IREF=(NUM(I,1)-1)*(IDIM+1)
              DO 310 J=1,I
                JREF=(NUM(J,1)-1)*(IDIM+1)
                VM=0.D0
                DO 320 I1=1,IDIM
                  PSCAL=0.D0
                  DO 330 J1=1,IDIM
                    PSCAL = PSCAL +
     *                       (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1)
  330             CONTINUE
                  VM=VM+(PSCAL*ZINVL(I1))**2
  320           CONTINUE
                VM=SQRT(VM)
                RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  310         CONTINUE
  300       CONTINUE
C
          ENDIF
C
        ENDIF
C
      ELSEIF (ILOI.EQ.2) THEN
C
C Loi gaussienne
C
        IF (IANI.EQ.1) THEN
C
C Cas isotrope
C
          DO 400 I=1,NBLI
            IREF=(NUM(I,1)-1)*(IDIM+1)
            DO 410 J=1,I
              JREF=(NUM(J,1)-1)*(IDIM+1)
              VM=0.D0
              DO 420 I1=1,IDIM
                VM=VM+(XCOOR(IREF+I1)-XCOOR(JREF+I1))**2
  420         CONTINUE
              RE(I,J,1)=ZSIG*ZSIG*EXP(-VM*ZINVL(1))
  410       CONTINUE
  400     CONTINUE
C
        ELSEIF (IANI.EQ.2) THEN
C
C Cas anisotrope :
C  directions d'anisotropie parallèles aux axes (IDIR nul)
C
          IF (IDIR.EQ.0) THEN
C
            DO 500 I=1,NBLI
              IREF=(NUM(I,1)-1)*(IDIM+1)
              DO 510 J=1,I
                JREF=(NUM(J,1)-1)*(IDIM+1)
                VM=0.D0
                DO 520 I1=1,IDIM
                  VM=VM+((XCOOR(IREF+I1)-XCOOR(JREF+I1))*ZINVL(I1))**2
  520           CONTINUE
                RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  510         CONTINUE
  500       CONTINUE
C
C Cas anisotrope :
C  directions d'anisotropie non parallèles aux axes (IDIR non nul)
C
          ELSEIF (IDIR.NE.0) THEN
C
            DO 600 I=1,NBLI
              IREF=(NUM(I,1)-1)*(IDIM+1)
              DO 610 J=1,I
                JREF=(NUM(J,1)-1)*(IDIM+1)
                VM=0.D0
                DO 620 I1=1,IDIM
                  PSCAL=0.D0
                  DO 630 J1=1,IDIM
                    PSCAL = PSCAL +
     *                       (XCOOR(IREF+J1)-XCOOR(JREF+J1))*XCOV(I1,J1)
  630             CONTINUE
                  VM=VM+(PSCAL*ZINVL(I1))**2
  620           CONTINUE
                RE(I,J,1)=ZSIG*ZSIG*EXP(-VM)
  610         CONTINUE
  600       CONTINUE
C
          ENDIF
C
        ENDIF
C
      ENDIF
      SEGDES,MCOORD
C
C Calcul de M telle que (RE =) COV = M Mt
C   M est stockée dans la partie inférieure de RE
C
      CALL SYCHOM(RE,NBLI)
C
C Création de l'objet rigidité
C
      CALL ACTOBJ('RIGIDITE',MRIGID,0)
      CALL ECROBJ('RIGIDITE',MRIGID)
C
      END