C BIOSAV SOURCE CB215821 25/03/25 21:15:04 12217 SUBROUTINE BIOSAV C---------------------------------------------------------- C C CALCUL DES CHAMPS DE BIOT ET SAVART C C---------------------------------------------------------- IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SMELEME -INC SMCOORD -INC TMTRAV -INC SMCHPOI -INC SMMODEL -INC SMCHAML INTEGER ENT DIMENSION C(3),PN1(3),RESP(3) DIMENSION VX(3),VZ(3),VY(3),VI(3),PP3(3) DIMENSION ZF1(1),ZALF(1),ZDENS(1),ZBET(1) DIMENSION P1(3),P2(3),P3(3),ROTA(3,3) SEGMENT ICPR(NBPTS) SEGMENT SANGLE REAL*8 TETM(NT4) REAL*8 TETI(NT4) REAL*8 DTE(NT4) ENDSEGMENT C CHARACTER*4 MCLE(4),CHSECT(1),CHOIX(2) CHARACTER*8 LISUP(5),MOT LOGICAL INDU,POTE,CHPO C DATA EPSIL3/1.D-12/ C DATA LISUP / 'NOEUD ','GRAVITE ','RIGIDITE','MASSE ', & 'STRESSES'/ DATA MCLE/'CERC','ARC ','BARR','FIL '/ DATA CHSECT/'TRAP'/ DATA CHOIX/'INDU','POTE'/ SEGACT MCOORD C SYNTAXE 2 : AVEC UN CHPOINT CALL LIROBJ('CHPOINT ',MCHPO1,0,IRETOU) IF (IERR.NE.0) RETURN IF (IRETOU.EQ.1) THEN CALL ACTOBJ('CHPOINT ',MCHPO1,1) IF (IERR.NE.0) RETURN CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU) IF (IERR.NE.0) RETURN CALL ACTOBJ('MAILLAGE',IPT1,1) IF (IERR.NE.0) RETURN C CALCUL PAR INTEGRALE ELLIPTIQUE (PHR 2002) CALL BIOKAM (MCHPO1,IPT1,ENT) IF (IERR.NE.0) RETURN MCHPOI=ENT C ATTRIBUTION D'UNE NATURE DISCRETE AU CHPOINT SEGACT MCHPOI*MOD JATTRI(1) = 2 C ECRITURE DANS LA PILE ET SORTIE CALL ACTOBJ('CHPOINT ',MCHPOI,1) CALL ECROBJ('CHPOINT ',MCHPOI) SEGDES MCOORD RETURN ENDIF C C SYNTAXE 1 : AVEC UN MAILLAGE OU UN MODELE C QUELQUES INITIALISATIONS GY = 0.D0 PENT1= 0.D0 PENT2= 0.D0 NB = 1 C CHOIX DE L'OPTION DE CALCUL (INDUCTION PAR DEFAUT) INDU=.TRUE. POTE=.FALSE. CALL LIRMOT(CHOIX,2,IMC,0) IF (IERR.NE.0) RETURN IF(IMC.EQ.2) THEN POTE=.TRUE. INDU=.FALSE. ENDIF C C LECTURE D'UN MAILLAGE OU BIEN D'UN MODELE CHPO=.FALSE. CALL LIROBJ('MAILLAGE',MELEME,0,IRETOU) IF (IERR.NE.0) RETURN IF (IRETOU.EQ.1) THEN CALL ACTOBJ('MAILLAGE',MELEME,1) IF (IERR.NE.0) RETURN CHPO=.TRUE. ELSE CALL LIROBJ('MMODEL',MMODE1,1,IRETOU) IF (IERR.NE.0) RETURN CALL ACTOBJ('MMODEL',MMODE1,1) IF (IERR.NE.0) RETURN C LECTURE FACULTATIVE D'UN LIEU SUPPORT CALL LIRMOT(LISUP,5,ISUP,0) IF (IERR.NE.0) RETURN C PAR DEFAUT, LE SUPPORT EST AUX NOEUDS IF (ISUP.EQ.0) ISUP=1 MOT=LISUP(ISUP) C ON CREE UN MCHAML NUL AUX POINTS SUPPORTS CALL ZEROP(MMODE1,MOT,MCHELM) C RECUPERATION DES CHAMPS DE COORDONNEES DES POINTS SUPPORTS: MCHEL1,2,3 CALL CHELCO(0,MCHELM,MCHEL3,MCHEL2,MCHEL1) IF (IDIM.EQ.2) THEN MCHEL1=MCHEL2 MCHEL2=MCHEL3 ENDIF C SUPPRESSION DU MCHAML NUL CALL DTCHAM(MCHELM) ENDIF C C LECTURE D'UN MOT CLE SUR LA FORME DE L'INDUCTEUR CALL LIRMOT(MCLE,4,IMC,1) IF (IERR.NE.0) RETURN C IMC = 1 ON A LU CERCLE C IMC = 2 ON A LU ARC C IMC = 3 ON A LU BARRE C IMC = 4 ON A LU FIL C C C SI 'CERC' EN 2D AXI, ON NE LIT QU'UN FLOTTANT IF ((IMC.EQ.1).AND.(IFOMOD.EQ.0)) THEN CALL LIRREE(Z1,1,IRETOU) C SINON, LECTURE DE 2 ou 3 POINTS ELSE CALL LIROBJ('POINT ',IP1,1,IRETOU) CALL LIROBJ('POINT ',IP2,1,IRETOU) IP3=IP2 C SI MOT CLEF AUTRE QUE 'FIL', ON LIT UN TROISIEME POINT IF(IMC.LT.4) THEN CALL LIROBJ('POINT ',IP3,1,IRETOU) C CAS DU FIL RECTILIGNE ELSE GOTO 10 ENDIF ENDIF IF (IERR.NE.0) RETURN C C LECTURE DE 4 OU 5 FLOTTANTS SELON LES CAS CALL LIRREE(F1,1,IRETOU) CALL LIRREE(F2,1,IRETOU) F3=0.D0 IF(IMC.LT.3) CALL LIRREE(F3,1,IRETOU) C C LECTURE DU TYPE DE SECTION ISECT=0 CALL LIRMOT(CHSECT,1,ISECT,0) IF(ISECT.EQ.1) THEN CALL LIRREE(PENT1,1,IRETOU) CALL LIRREE(PENT2,1,IRETOU) IF(IERR.NE.0) RETURN ELSE ISECT=3 ENDIF 10 CONTINUE C LECTURE DE LA DENSITE DE COURANT ET DU LA PERMEABILITE CALL LIRREE(DENS,1,IRETOU) CALL LIRREE(XIU,1,IRETOU) C DANS LE CAS 'CERC' EN 2D AXI, LA DENSITE DE COURANT ORIENTEE C SELON UT = UR ^ UZ IF ((IMC.EQ.1).AND.(IFOMOD.EQ.0)) THEN DENS=-DENS ENDIF C ON VIENT DE LIRE MU0 DANS L UNITE UTILISATEUR RUMKSA=4.D0*XPI*1.D-7 RAP=XIU/RUMKSA IU=1 IF (IERR.NE.0) RETURN ZF1(1)=F1 ZDENS(1)=DENS C RECUPERATION DES COORDONNEES DU POINT P1 C SI 2D AXI, ON PREND P1=(X=0,Y=0,Z=Z1) IF (IFOMOD.EQ.0) THEN P1(1)=0.D0 P1(2)=0.D0 P1(3)=Z1 ELSE IREF=(IDIM+1)*IP1 P1(1)=XCOOR(IREF-IDIM) P1(2)=XCOOR(IREF-IDIM+1) P1(3)=0.D0 IF (IDIM.EQ.3) P1(3)=XCOOR(IREF-IDIM+2) ENDIF C RECUPERATION DES COORDONNEES DU POINT P2 C SI 2D AXI, ON PREND P2=(X=F1,Y=0,Z=Z1) IF (IFOMOD.EQ.0) THEN P2(1)=F1 P2(2)=0.D0 P2(3)=Z1 C SINON, ON RECUPERE LES COORDONNEES DU SECOND POINT LU ELSE IREF=(IDIM+1)*IP2 P2(1)=XCOOR(IREF-IDIM) P2(2)=XCOOR(IREF-IDIM+1) P2(3)=0.D0 IF (IDIM.EQ.3) P2(3)=XCOOR(IREF-IDIM+2) ENDIF C RECUPERATION DES COORDONNEES DU POINT P3 (SI MOT CLEF AUTRE QUE 'FIL') IF(IMC.LT.4) THEN C SI 2D AXI, ON PREND P3=(X=0,Y=F1,Z=Z1) IF (IFOMOD.EQ.0) THEN P3(1)=0.D0 P3(2)=F1 P3(3)=Z1 C SINON, ON RECUPERE LES COORDONNEES DU TROISIEME POINT LU ELSE IREF=(IDIM+1)*IP3 P3(1)=XCOOR(IREF-IDIM) P3(2)=XCOOR(IREF-IDIM+1) P3(3)=0.D0 IF (IDIM.EQ.3) P3(3)=XCOOR(IREF-IDIM+2) ENDIF ENDIF C C ----------------------- C | LES CALCULS .... | C ----------------------- C C ON COMMENCE PAR EXPLORER LE MAILLAGE OU LE MODELE C ET INITIALISER LE CHAMP RESULTAT C C CAS DES MAILLAGES (RESULTAT CHPOINT) IF (CHPO) THEN SEGACT MELEME C CONVERTION EN 'POI1' SI NECESSAIRE IF (ITYPEL.NE.1.OR.LISOUS(/1).NE.0) CALL CHANGE(MELEME,1) C ON INITIALISE LE SEGMENT DE TRAVAIL DU CHPOINT RESULTAT C NOMS DES HARMONIQUES ET DES COMPOSANTES C POUR L'INSTANT LIMITE AU 3D, 2D PLAN ET 2D AXIS NNNOE=NUM(/2) NNIN=IDIM SEGINI MTRAV NHRM=NIFOUR IF(NIFOUR.EQ.1) NHRM=IFOUR NHAR(1)=NHRM NHAR(2)=NHRM IF (IDIM.EQ.3) NHAR(3)=NHRM C SI 2D AXI, LES COMPOSANTES SERONT SELON R,Z IF (IFOMOD.EQ.0) THEN IF (INDU) THEN INCO(1)='BR ' INCO(2)='BZ ' ELSE INCO(1)='AR ' INCO(2)='AZ ' ENDIF C SINON, LES COMPOSANTES SONT SELON X,Y,Z ELSE IF (INDU) THEN INCO(1)='BX ' INCO(2)='BY ' INCO(3)='BZ ' ELSE INCO(1)='AX ' INCO(2)='AY ' INCO(3)='AZ ' ENDIF ENDIF GOTO (100,100,200,300),IMC C C CAS DU MODELE (RESULTAT MCHAML) ELSE C PREPARATION DU CHAMP RESULTAT SUR LA BASE DE CELUI DE LA COORDONNEE 1 : MCHEL1 N1=MCHEL1.INFCHE(/1) N3=MCHEL1.INFCHE(/2) L1=9 SEGADJ MCHEL1 IF(INDU) THEN MCHEL1.TITCHE='INDUCTION' ELSE MCHEL1.TITCHE='POTENTIEL' ENDIF C BOUCLE SUR LES SOUS ZONES NSZ=N1 DO ISZ=1,NSZ MCHAM1=MCHEL1.ICHAML(ISZ) MCHAM2=MCHEL2.ICHAML(ISZ) IF (IDIM.EQ.3) MCHAM3=MCHEL3.ICHAML(ISZ) C AJOUT DE COMPOSANTES DANS MCHAM1 N2=IDIM SEGADJ MCHAM1 C RENOMAGE DES COMPOSANTES C SI 2D AXI, LES COMPOSANTES SERONT SELON R,Z IF (IFOMOD.EQ.0) THEN IF (INDU) THEN MCHAM1.NOMCHE(1)='BR ' MCHAM1.NOMCHE(2)='BZ ' ELSE MCHAM1.NOMCHE(1)='AR ' MCHAM1.NOMCHE(2)='AZ ' ENDIF C SINON, LES COMPOSANTES SONT SELON X,Y,Z ELSE IF (INDU) THEN MCHAM1.NOMCHE(1)='BX ' MCHAM1.NOMCHE(2)='BY ' MCHAM1.NOMCHE(3)='BZ ' ELSE MCHAM1.NOMCHE(1)='AX ' MCHAM1.NOMCHE(2)='AY ' MCHAM1.NOMCHE(3)='AZ ' ENDIF ENDIF C AJOUT DES MELVAL DES AUTRES COORDONNEES MCHAM1.TYPCHE(2)=MCHAM2.TYPCHE(1) MCHAM1.IELVAL(2)=MCHAM2.IELVAL(1) SEGSUP MCHAM2 IF (IDIM.EQ.3) THEN MCHAM1.TYPCHE(3)=MCHAM3.TYPCHE(1) MCHAM1.IELVAL(3)=MCHAM3.IELVAL(1) SEGSUP MCHAM3 ENDIF ENDDO SEGSUP MCHEL2 IF (IDIM.EQ.3) SEGSUP MCHEL3 GOTO (100,100,200,300),IMC ENDIF 100 CONTINUE C BOBINE CIRCULAIRE C ETABLISSEMENT DU REPERE LOCAL C AXE DES Z (P2 - P1 ) ^ (P3 - P1) DO I=1,3 VX(I)=P2(I)-P1(I) VI(I)=P3(I)-P1(I) ENDDO CALL PROVEC(VX,VI,VZ) CALL PROVEC(VZ,VX,VY) CALL ROTREP (ROTA,VX,VY,VZ) IF(IERR.NE.0) THEN IF (CHPO) SEGSUP MTRAV RETURN ENDIF IF(A_EGALE_B(F1,XZERO)) THEN SEGDES MCOORD RETURN ENDIF ALF=F2/F1 BET=F3/(2.D0*F1) ZALF(1)=ALF ZBET(1)=BET C IC=3 RAYO2=EPSIL3 PROD=0.D0 DO J=1,3 PP3(J)= P3(J)-P1(J) RAYO2 = RAYO2+VX(J)*VX(J) PROD = PROD+VX(J)*VI(J) ENDDO TM=ACOS(PROD/RAYO2) GOTO 350 200 CONTINUE C CAS DES BARRES ON VA TOURNER / OX POUR AMENER OY // P1P2 DO I=1,3 VX(I)=P2(I)-P1(I) VI(I)=P3(I)-P1(I) GY=GY+(VX(I)*VX(I)) C(I)= (P2(I)+ P1(I))*0.5D0 ENDDO GY=SQRT( GY) CALL PROVEC(VX,VI,VZ) CALL PROVEC(VZ,VX,VY) CALL ROTREP (ROTA,VX,VY,VZ) IF(IERR.NE.0) THEN IF (CHPO) SEGSUP MTRAV RETURN ENDIF DO I= 1,3 VZ(I)=0.D0 DO J=1,3 VZ(I)=VZ(I)+ROTA(I,J)*C(J) ENDDO ENDDO GOTO 350 C CAS DU FIL 300 CONTINUE 350 CONTINUE UPSI =1.D-9 C INITIALISATION DU SEGMENT SANGLE NT4=8000 SEGINI SANGLE C CAS AVEC MAILLAGE : ON BOUCLE SUR LES NOEUDS IF(CHPO) THEN DO IPT=1,NNNOE C L EST LE NUMERO DU NOEUD DANS LE TABLEAU DES COORDONNEES C RECUPERATION DES COORDONNEES DU NOEUD OU L'ON VA CALCULER LE CHAMP C DANS LE TABLEAU PN1 L=NUM(1,IPT) IREF=(IDIM+1)*L DO J=1,3 PN1(J)=XCOOR(IREF-IDIM-1+J) ENDDO C SI 2D AXI, ON PREND PN1=(X=R,Y=0,Z=Z) IF (IFOMOD.EQ.0) THEN PN1(3)=PN1(2) PN1(2)=0.D0 ENDIF C CALCUL DU CHAMP (B OU A) AU POINT PN1 CALL BIOSA1(INDU,POTE,IMC,PN1,P1,P2,PP3,ROTA,VI,GY,IU,NB,TM, & UPSI,PENT1,PENT2,F1,F2,F3,ZF1,ZALF,ZDENS,ZBET,C, & ISECT,SANGLE,VX,VY,VZ,RESP) IF(IERR.NE.0) THEN SEGSUP MTRAV RETURN ENDIF C RETOUR DU CHAMP DANS LE REPERE GENERAL DO I=1,3 VY(I)=0.D0 DO J=1,3 VY(I)=VY(I)+ROTA(J,I)*RESP(J)*RAP ENDDO ENDDO C SI 2D AXI, ON VA CHERCHER LA COMPOSANTE 3 ET ON LA MET EN POSITION 2 C POUR AVOIR : BR,BZ (OU AR,AZ) IF (IFOMOD.EQ.0) VY(2)=VY(3) C REMPLISSAGE DU SEGMENT DE TRAVAIL MTRAV C REFERENCE DU NOEUD, VALEUR DU CHAMP ET INDICATEUR DE DEFINITION IGEO(IPT)=L DO J=1,IDIM BB(J,IPT)=VY(J) IBIN(J,IPT)=1 ENDDO ENDDO C CAS AVEC MMODELE : ON BOUCLE SUR LES SOUS ZONES/ELEMENTS/POINTS SUPPORTS ELSE DO ISZ=1,NSZ C RECUPERATION DES COORDONNEES DU POINT SUPPORT OU L'ON VA CALCULER LE CHAMP C DANS LE TABLEAU PN1 MCHAM1=MCHEL1.ICHAML(ISZ) MELVA1=MCHAM1.IELVAL(1) MELVA2=MCHAM1.IELVAL(2) IF (IDIM.EQ.3) MELVA3=MCHAM1.IELVAL(3) NP=MELVA1.VELCHE(/1) NE=MELVA1.VELCHE(/2) DO IE=1,NE DO K=1,NP PN1(1)=MELVA1.VELCHE(K,IE) PN1(2)=MELVA2.VELCHE(K,IE) PN1(3)=0.D0 IF (IDIM.EQ.3) PN1(3)=MELVA3.VELCHE(K,IE) C SI 2D AXI, ON PREND PN1=(X=R,Y=0,Z=Z) IF (IFOMOD.EQ.0) THEN PN1(3)=PN1(2) PN1(2)=0.D0 ENDIF C CALCUL DU CHAMP (B OU A) AU POINT PN1 CALL BIOSA1(INDU,POTE,IMC,PN1,P1,P2,PP3,ROTA,VI,GY,IU,NB, & TM,UPSI,PENT1,PENT2,F1,F2,F3,ZF1,ZALF,ZDENS, & ZBET,C,ISECT,SANGLE,VX,VY,VZ,RESP) IF(IERR.NE.0) THEN SEGSUP MCHEL1 RETURN ENDIF C RETOUR DU CHAMP DANS LE REPERE GENERAL DO II=1,3 VY(II)=0.D0 DO JJ=1,3 VY(II)=VY(II)+ROTA(JJ,II)*RESP(JJ)*RAP ENDDO ENDDO C SI 2D AXI, ON VA CHERCHER LA COMPOSANTE 3 ET ON LA MET EN POSITION 2 C POUR AVOIR : BR,BZ (OU AR,AZ) IF (IFOMOD.EQ.0) VY(2)=VY(3) C ON REMPLACE LES VALEURS DES COORDONNEES DU POINT PAR CELLE DU CHAMP CALCULE MELVA1.VELCHE(K,IE)=VY(1) MELVA2.VELCHE(K,IE)=VY(2) IF (IDIM.EQ.3) MELVA3.VELCHE(K,IE)=VY(3) ENDDO ENDDO ENDDO ENDIF C SUPRESSION DU SEGMENT SANGLE SEGSUP SANGLE C CREATION DU CHAMP RESULTAT C SI CHPOINT, CELUI-CI EST FAIT A PARTIR DU SEGMENT DE TRAVAIL MTRAV IF (CHPO) THEN CALL CRECHP(MTRAV,MCHPOI) SEGSUP MTRAV C ATTRIBUTION D'UNE NATURE DISCRETE AU CHPOINT SEGACT MCHPOI*MOD JATTRI(1)=2 C ECRITURE DANS LA PILE CALL ACTOBJ('CHPOINT ',MCHPOI,1) CALL ECROBJ('CHPOINT ',MCHPOI) C SI MCHAML ELSE CALL ACTOBJ('MCHAML ',MCHEL1,1) CALL ECROBJ('MCHAML ',MCHEL1) ENDIF SEGDES MCOORD END