C TGRAD SOURCE CHAT 09/10/09 21:25:20 6519 SUBROUTINE TGRAD (IPTAB1) C C*********************************************************************** C C C FONCTION: C --------- C C en "0D/1D" sur des éléments de type POINT. C C C ENTREE : C -------- C C IPTAB1 : TABLE de soustype 'OPER_0D' contenant les indices C suivants : C C TAB1 . 'GEOINF' : TABLE des informations géométriques de soustype C 'GEOINF' (type ENTIER). C TAB1 . 'INCO' : TABLE de soustype 'INCO' contenant l'ensemble C des champs à l'itération précédant l'itération C courante (type ENTIER). C TAB1 . 'POTENTIA' : Valeurs du champ potentiel (gravité * cote) C (type CHPOINT) (indice facultatif). C Cet indice est utilisé pour calculer p.e le C gradient de rho * gravité * cote. C (doit avoir le même support que l'inconnue C duale). C TAB1 . 'DUAL' : Nom de l'inconnue duale (doit être un indice de C la table de soustype 'INCO') (type MOT). C TAB1 . 'PRIMAL' : Nom de l'inconnue primale (doit être un indice C de la table de soustype 'INCO') (type MOT). C C C RESULTAT : C ---------- C C TAB1 . 'LHS' : Matrice élémentaire associée à l'opération C (type RIGIDITE). C C C AUTEUR, DATE DE CREATION: C ------------------------- C C Laurent DADA décembre 1996 C C C LANGAGE: C -------- C C ESOPE + FORTRAN77 C C*********************************************************************** C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C -INC CCOPTIO -INC CCGEOME -INC SMCOORD -INC SMTABLE POINTEUR MTABG.MTABLE,MTABS.MTABLE,IPTABI.MTABLE -INC SMCHPOI -INC SMELEME POINTEUR IPTD1.MELEME,IPTP1.MELEME,IPJUCE.MELEME -INC SMRIGID C SEGMENT REDI INTEGER ORDR1(NNGOT) INTEGER ORDR2(NNGOT) ENDSEGMENT C CHARACTER*8 TYPE,MOTI,MOT1,NOMPR1,NOMDU1,NOSUP1,NOSUD1 CHARACTER*7 NAMT1 LOGICAL LPOT C C Lecture de la table GEOINF de la table OPER_0D C TYPE = 'TABLE ' CALL ACMO (IPTAB1,'GEOINF',TYPE,IPTABG) IF (IERR.NE.0) RETURN C MOTI = 'SOUSTYPE' CALL ACMM (IPTABG,MOTI,MOT1) IF (IERR.NE.0) RETURN IF (MOT1(1:6).NE.'GEOINF') THEN MOTERR(1:8) = 'GEOINF ' MOTERR(9:16) = 'GEOINF ' CALL ERREUR (790) RETURN ENDIF C C Lecture de la table INCO dans la table OPER_0D C TYPE = 'TABLE ' CALL ACMO (IPTAB1,'INCO',TYPE,IPTABI) IF (IERR.NE.0) RETURN C MOTI = 'SOUSTYPE' CALL ACMM (IPTABI,MOTI,MOT1) IF (IERR.NE.0) RETURN IF (MOT1(1:4).NE.'INCO') THEN MOTERR(1:8) = 'INCO ' MOTERR(9:16) = 'INCO ' CALL ERREUR (790) RETURN ENDIF C C Lecture de la table SUPPORT dans la table INCO C TYPE = 'TABLE ' CALL ACMO (IPTABI,'SUPPORT',TYPE,IPTABS) IF (IERR.NE.0) RETURN C C Lecture du MAILLAGE des connectivités 'JUNCEL' de la table GEOINF C Arrêt si les éléments ne sont pas des SEG3 C TYPE = 'MAILLAGE' CALL ACMO (IPTABG,'JUNCEL',TYPE,IPJUCE) IF (IERR.NE.0) RETURN SEGACT IPJUCE IF ((IPJUCE.ITYPEL).NE.3) THEN MOTERR(1:8) = 'JUNCEL ' MOTERR(9:16) = 'MAILLAGE' CALL ERREUR (787) SEGDES IPJUCE RETURN ENDIF NBEJC1 = IPJUCE.NUM(/2) C C Lecture du nom de l'inconnue PRIMAL C TYPE = ' ' CALL ACMO (IPTAB1,'PRIMAL',TYPE,IPR1) IF (IERR.NE.0) RETURN IF (TYPE.EQ.'MOT ') THEN CALL ACMM (IPTAB1,'PRIMAL',NOMPR1) IF (IERR.NE.0) RETURN ELSE MOTERR(1:8) = 'PRIMAL ' MOTERR(9:16) = TYPE CALL ERREUR (787) RETURN ENDIF C C Lecture du nom de l'inconnue DUAL C TYPE = ' ' CALL ACMO (IPTAB1,'DUAL',TYPE,IDU1) IF (IERR.NE.0) RETURN IF (TYPE.EQ.'MOT ') THEN CALL ACMM (IPTAB1,'DUAL',NOMDU1) IF (IERR.NE.0) RETURN ELSE MOTERR(1:8) = 'DUAL ' MOTERR(9:16) = TYPE CALL ERREUR (787) RETURN ENDIF C C Lecture éventuelle du champ POTENTIA C LPOT = .FALSE. TYPE = ' ' CALL ACMO (IPTAB1,'POTENTIA',TYPE,IPCHP1) IF (IERR.NE.0) RETURN IF (TYPE.EQ.'CHPOINT ') THEN LPOT = .TRUE. C Contrôle du CHPOINT MCHPOI = IPCHP1 SEGACT MCHPOI MSOUPO = IPCHP(1) SEGDES MCHPOI SEGACT MSOUPO NC = NOHARM(/1) IF (NC.NE.1) THEN MOTERR(1:8) = 'POTENTIA' MOTERR(9:16) = 'CHPOINT ' CALL ERREUR (784) SEGDES MSOUPO RETURN ENDIF IPT1 = IGEOC MPOVA1 = IPOVAL SEGDES MSOUPO ENDIF C C Lecture du nom du support de l'inconnue PRIMAL C Lecture du MAILLAGE de l'inconnue PRIMAL C TYPE = ' ' CALL ACMO (IPTABS,NOMPR1,TYPE,ISUP1) IF (IERR.NE.0) RETURN IF (TYPE.EQ.'MOT ') THEN CALL ACMM (IPTABS,NOMPR1,NOSUP1) IF (IERR.NE.0) RETURN ENDIF TYPE = 'MAILLAGE' CALL ACMO (IPTABG,NOSUP1,TYPE,IPTP1) IF (IERR.NE.0) RETURN C C Lecture du nom du support de l'inconnue DUAL C Arrêt si différent de 'JUNCTION' C Lecture du MAILLAGE de l'inconnue DUAL C TYPE = ' ' CALL ACMO (IPTABS,NOMDU1,TYPE,ISUD1) IF (TYPE.EQ.'MOT ') THEN CALL ACMM (IPTABS,NOMDU1,NOSUD1) IF (IERR.NE.0) RETURN IF (NOSUD1.NE.'JUNCTION') THEN MOTERR(1:8) = 'DUAL ' MOTERR(9:16) = 'CHPOINT ' CALL ERREUR (788) RETURN ENDIF ENDIF TYPE = 'MAILLAGE' CALL ACMO (IPTABG,NOSUD1,TYPE,IPTD1) IF (IERR.NE.0) RETURN C C Création du support géométrique pour la RIGIDITE C (maillage de type SUPER-ELEMENT). C NNGOT = XCOOR(/1)/(IDIM+1) SEGINI REDI C C On fusionne les maillages de POI1 des supports des inconnues PRIMAL C et DUAL en un maillage de type SUPER-ELEMENT. C SEGACT IPTD1 IF (IPTD1.ITYPEL.NE.1) CALL CHANGE (IPTD1,1) SEGACT IPTD1 SEGACT IPTP1 IF (IPTP1.ITYPEL.NE.1) CALL CHANGE (IPTP1,1) SEGACT IPTP1 NBNNP1 = IPTP1.NUM(/2) NBNND1 = IPTD1.NUM(/2) NBNN = NBNNP1 + NBNND1 NBSOUS = 0 NBREF = 0 NBELEM = 1 SEGINI MELEME ICOLOR(1) = IDCOUL ITYPEL = 28 DO 50 I50=1,NBNNP1 NUM(I50,1) = IPTP1.NUM(1,I50) ORDR2(NUM(I50,1)) = I50 50 CONTINUE DO 60 I60=1,NBNND1 NUM(I60+NBNNP1,1) = IPTD1.NUM(1,I60) ORDR1(NUM(I60+NBNNP1,1)) = I60 60 CONTINUE C C Création de la RIGIDITE C NRIGE = 7 NRIGEL = 1 SEGINI MRIGID C MTYMAT = 'RIGIDITE' IFORIG = IFOMOD ICHOLE = 0 IMGEO1 = 0 IMGEO2 = 0 ISUPEQ = 0 COERIG(1) = 1.D0 IRIGEL(1,1) = MELEME IRIGEL(2,1) = 0 IRIGEL(5,1) = NIFOUR IRIGEL(6,1) = 0 IRIGEL(7,1) = 2 C SEGDES MELEME C C Remplissage du descripteur de l'objet RIGIDITE C NLIGRP = NBNNP1 NLIGRD = NBNND1 SEGINI DESCR IRIGEL(3,1) = DESCR DO 10 I10=1,NBNNP1 NOELEP(I10) = I10 LISINC(I10) = NOMPR1 10 CONTINUE DO 11 I11=1,NBNND1 NOELED(I11) = I11+NBNNP1 LISDUA(I11) = NOMDU1 11 CONTINUE C SEGDES DESCR C NELRIG = 1 * SEGINI IMATRI SEGINI XMATRI C * IMATTT(1) = XMATRI IRIGEL(4,1) = xMATRI * SEGDES IMATRI C C Remplissage de la matrice élémentaire C C mise à 0.d0 de la matrice élémentaire DO 70 I70=1,NBNND1 DO 701 I701=1,NBNNP1 RE(I70,I701,1) = 0.D0 701 CONTINUE 70 CONTINUE C C avec prise en compte du terme POTENTIA C IF (LPOT) THEN C activation du maillage et des valeurs du CHPOINT POTENTIA SEGACT IPT1 SEGACT MPOVA1 C petit controle du support du champ POTENTIA NBEL1 = IPT1.NUM(/2) IF (NBEL1.NE.(NBNNP1+NBNND1)) THEN MOTERR(1:8) = 'POTENTIA' MOTERR(9:16) = 'CHPOINT ' CALL ERREUR (788) SEGDES IPT1 SEGDES MPOVA1 SEGDES IPJUCE SEGSUP MELEME SEGSUP XMATRI SEGSUP MRIGID SEGDES IPTP1 SEGDES IPTD1 SEGSUP REDI RETURN ENDIF DO 72 I72=1,NBEJC1 NPT1 = IPJUCE.NUM(1,I72) NPTF1 = IPJUCE.NUM(2,I72) NPT2 = IPJUCE.NUM(3,I72) C récupération des valeurs du champ POTENTIA aux points NPT1,NPT2,NPTF1 C balayage sur les points du CHPOINT POTENTIA XGZP1 = 0.D0 XGZD1 = 0.D0 XGZP2 = 0.D0 DO 723 I723=1,NBEL1 IF (NPT1.EQ.(IPT1.NUM(1,I723))) XGZP1 = MPOVA1.VPOCHA(I723,1) IF (NPT2.EQ.(IPT1.NUM(1,I723))) XGZP2 = MPOVA1.VPOCHA(I723,1) IF (NPTF1.EQ.(IPT1.NUM(1,I723))) XGZD1 = MPOVA1.VPOCHA(I723,1) 723 CONTINUE C positions dans la matrice élémentaire IPRI1 = ORDR2(NPT1) IDUA1 = ORDR1(NPTF1) IPRI2 = ORDR2(NPT2) C remplissage de la matrice RE(IDUA1,IPRI1,1) = -1.D0 * (XGZP1 - XGZD1) RE(IDUA1,IPRI2,1) = (XGZP2 - XGZD1) 72 CONTINUE SEGDES IPTP1 SEGDES IPTD1 C C ELSE C C DO 82 I82=1,NBEJC1 NPT1 = IPJUCE.NUM(1,I82) NPTF1 = IPJUCE.NUM(2,I82) NPT2 = IPJUCE.NUM(3,I82) C positions dans la matrice élémentaire IPRI1 = ORDR2(NPT1) IDUA1 = ORDR1(NPTF1) IPRI2 = ORDR2(NPT2) C remplissage de la matrice RE(IDUA1,IPRI1,1) = -1.D0 RE(IDUA1,IPRI2,1) = 1.D0 82 CONTINUE ENDIF C SEGDES IPJUCE SEGDES XMATRI SEGSUP REDI C C Ecriture du résultat C TYPE = 'RIGIDITE' CALL ECMO (IPTAB1,'LHS',TYPE,MRIGID) IF (IERR.NE.0) RETURN C SEGDES MRIGID C END