reduaf
C REDUAF SOURCE PV090527 24/04/22 21:15:04 11913 C Reduction du champ par element jchelm sur le modele mmodtm C Le resultat est le champ par element mchel2 pour iret = 1 (KERRE=0), C sinon en cas d'erreur mchel2 = 0 pour iret = 0 (KERRE = num. erreur) C En sortie le champ mchel2 est un champ entierement actif. C (kich) en sortie mmodel deroule IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER (I-N) -INC PPARAM -INC CCOPTIO -INC SMCHAML -INC SMMODEL -INC SMCOORD -INC SMELEME -INC SMLENTI -INC CCPRECO -INC CCASSIS EXTERNAL LONG segment izone(NZ,NSMOD) segment ismel(NZ,NSMOD) segment szsxx logical lzsxx(NZ,NSMOD) endsegment segment icpr(nbpt) segment inde(jg) CHARACTER*(NCONCH) conloc,MO24 CHARACTER*(LOCOMP) nomloc CHARACTER*(16) typloc CHARACTER*(50) typ1 LOGICAL BVALID,OOOVP1,dmopa melpv = 0 ith1 = oothrd + 1 CG if (iimpi.eq.7203) then CG write(ioimp,*) 'Entree dans reduaf',mmodtm,jchelm CG call zpchel(jchelm,1) CG endif c write(6,*) 'Entree dans reduaf',jchelm,istri,mmodtm CALL oooho2(ihcour) iret = 1 KERRE = 0 mchel2 = 0 MO24 =' ' IF(istri .EQ. 0)THEN C Extension du MMODEL en cas de modele de MELANGE et NON-STRICT ELSEIF(istri .EQ. 1)THEN mmodel = mmodtm ELSE ENDIF NSMOD = mmodel.kmodel(/1) DO is = 1, NSMOD imodel = mmodel.kmodel(is) C Verification si on a un modele de DARCY actuellement incompatible C Car il se servent du MAILLAGE dans la TABLE DOMAINE et pas celui C contenu dans le MMODEL IF (IDARC .NE. 0) THEN mchel2=jchelm RETURN ENDIF ENDDO C --------------------------------------------------------------------- C Preconditionnement de REDU C Verification que le resultat n'est pas deja dans le CCPRECO C --------------------------------------------------------------------- ITAILL = NBPRRE(ith1) CALL oooho1(mmodel,IHOmmo) CALL oooho1(jchelm,IHOjch) DO 201 IPREC1 = 1, ITAILL IF (PRECMO(IPREC1,ith1) .NE. mmodel) GOTO 201 IF ((PRECM1(IPREC1,ith1) .NE. jchelm) .AND. & (PRECM2(IPREC1,ith1) .NE. jchelm)) GOTO 201 IF (PRECM3(IPREC1,ith1) .NE. istri ) GOTO 201 C Ajout test horodatage du MMODEL et MCHAML d'entree (il a pu etre supprime puis recree avec le meme descripteur) IF (PRECM4(IPREC1,ith1) .NE. IHOmmo) GOTO 201 IF (PRECM5(IPREC1,ith1) .NE. IHOjch) GOTO 201 mchel2 = PRECM2(IPREC1,ith1) C IF (IPREC1 .EQ. NPREDU) THEN C PRINT *,' CCPRECO trop petit :',IPREC1 C CALL ERREUR(5) C ENDIF C Mise a jour du preconditionnement dans CCPRECO : Deplacement en position 1 du REDU deja fait IF (IPREC1 .EQ. 1) THEN RETURN ELSE DO IPREC2 = IPREC1,2,-1 PRECMO(IPREC2,ith1) = PRECMO(IPREC2 - 1,ith1) PRECM1(IPREC2,ith1) = PRECM1(IPREC2 - 1,ith1) PRECM2(IPREC2,ith1) = PRECM2(IPREC2 - 1,ith1) PRECM3(IPREC2,ith1) = PRECM3(IPREC2 - 1,ith1) PRECM4(IPREC2,ith1) = PRECM4(IPREC2 - 1,ith1) PRECM5(IPREC2,ith1) = PRECM5(IPREC2 - 1,ith1) ENDDO PRECMO(1,ith1) = mmodel PRECM1(1,ith1) = jchelm PRECM2(1,ith1) = mchel2 PRECM3(1,ith1) = istri PRECM4(1,ith1) = IHOmmo PRECM5(1,ith1) = IHOjch RETURN ENDIF 201 CONTINUE C 1 CONTINUE C Mise a jour du preconditionnement dans CCPRECO : Deplacement pour ecrire le nouveau REDU en 1ere position ITAILL = MIN(ITAILL + 1, NPREDU) NBPRRE(ith1) = ITAILL DO IPRECO = ITAILL,2,-1 PRECMO(IPRECO,ith1) = PRECMO(IPRECO - 1,ith1) PRECM1(IPRECO,ith1) = PRECM1(IPRECO - 1,ith1) PRECM2(IPRECO,ith1) = PRECM2(IPRECO - 1,ith1) PRECM3(IPRECO,ith1) = PRECM3(IPRECO - 1,ith1) PRECM4(IPRECO,ith1) = PRECM4(IPRECO - 1,ith1) PRECM5(IPRECO,ith1) = PRECM5(IPRECO - 1,ith1) ENDDO PRECMO(1,ith1) = mmodel PRECM1(1,ith1) = jchelm C PRECM2 doit etre mis a jour plus loin avant chaque RETURN PRECM3(1,ith1) = istri PRECM4(1,ith1) = IHOmmo PRECM5(1,ith1) = IHOjch mchelm = jchelm NZ = mchelm.imache(/1) L1 = mchelm.titche(/1) N3 = mchelm.infche(/2) C ----------------------------------------- C Cas tres particulier de MCHELM resultat : C ----------------------------------------- IF (NZ.EQ.0) THEN CG if (iimpi.eq.7203) write(ioimp,*) 'CAS PARTICULIER NZ = 0' C Mise a jour du preconditionnement dans CCPRECO mchel2 = jchelm PRECM2(1,ith1) = jchelm RETURN ENDIF C Quelques initialisations : C mlent2 contient le nombre d'elements du maillage de chaque sous-modele. jg = NSMOD call oooprl(1) SEGINI,mlent2,izone,ismel,szsxx C mlent3 contient les intersections entre les maillages determinees : C mlent3.lect(i3) avec ismel(iz,is) = i3 correspond a l'intersection C entre le maillage du sous-modele is et la sous-zone iz du champ si C la valeur de i3 n'est pas nulle ! jg = NSMOD * NZ SEGINI,mlent3 call oooprl(0) NL3 = 0 ISOZM = 0 icpr = 0 inde = 0 C C Regroupement des zones directement appariees avec un sous-modele C Recherche des zones pouvant intersecter le maillage d'un sous-modele DO 100 is = 1, NSMOD imodel = mmodel.kmodel(is) IF (imodel.nefmod.EQ.259) GOTO 100 meleme = imodel.imamod CALL oooho1(meleme,IHO1) itypm = meleme.itypel mlent2.lect(is) = meleme.num(/2) C On parcourt tous les NZ chamelem elementaires. DO 101 iz = 1, NZ conloc = mchelm.conche(iz) lzsxx(iz,is) = .false. IF (conloc.NE.MO24 .AND. & conloc .NE.imodel.conmod(1:LCONMO)) GOTO 101 C PRINT *,'REDUAF:',is,iz,':',conloc,':',imodel.conmod,':' ixx = 0 ipt1 = mchelm.imache(iz) C Correspondance maillage sous-zone et sous-modele IF (ipt1.EQ.meleme) THEN ixx = 1 lzsxx(iz,is) = .true. C Pas de correspondance directe, recherche intersection potentielle ELSE IF (ipt1.itypel.NE.itypm) GOTO 102 CALL oooho1(ipt1,IHO2) C Verification dans le PRECONDITIONNEMENT si deja evaluee DO 400 III=1,NINTSA(ith1) IF(PMAMOD(III,ith1) .NE. meleme) GOTO 400 IF(PMAMOH(III,ith1) .NE. IHO1 ) GOTO 400 IF(PMACHA(III,ith1) .NE. ipt1 ) GOTO 400 IF(PMACHH(III,ith1) .NE. IHO2 ) GOTO 400 mlenti=PMLENT(III,ith1) C PRINT *,'REDUAF_PRECONDITION',oothrd,meleme,ipt1,mlenti C IF(mlenti .EQ. 0) THEN C ixx = 0 C ismel(iz,is) = 0 C C ELSE NL3 = NL3 + 1 mlent3.lect(NL3) = mlenti ixx = -1 ismel(iz,is) = NL3 C ENDIF GOTO 102 400 CONTINUE C PRINT *,'REDUAF_INTERSECTION',oothrd,meleme,ipt1 C On va regarder si on n a pas deja evalue l'intersection : C (meme sous-modele is et sous-zone precedente ia<iz) DO ia = 1, iz-1 IF (ipt1.EQ.mchelm.imache(ia)) THEN IF (ismel(ia,is).GT.0) THEN ixx = -2 ismel(iz,is) = ismel(ia,is) GOTO 102 ENDIF ENDIF ENDDO C (meme sous-zone iz et sous-modele ia<is) DO 103 ia = 1, is-1 imode2 = mmodel.kmodel(ia) IF (imode2.nefmod.EQ.259) GOTO 103 ipt2 = imode2.imamod IF (ipt2.EQ.meleme) THEN IF (ismel(iz,ia).GT.0) THEN ixx = -3 ismel(iz,is) = ismel(iz,ia) GOTO 102 ENDIF ENDIF 103 CONTINUE C Détermination de l'intersection de ipt1 et meleme : C Creation d'un tableau (LISTENTI) de correspondance des C elements de IPT1 qui sont dans MELEME nbno1 = ipt1.num(/1) nbel1 = ipt1.num(/2) IF (icpr.EQ.0) THEN nbpt = nbpts + 1 np1 = nbpt - 1 SEGINI,icpr ELSE DO j = 1, nbpt icpr(j) = 0 ENDDO ENDIF DO j = 1, nbel1 DO m = 1, nbno1 ib = ipt1.num(m,j) icpr(ib) = icpr(ib) + 1 ENDDO ENDDO iprec = icpr(1) DO j = 2, np1 iprec = iprec + icpr(j) icpr(j) = iprec ENDDO jg = icpr(np1) icpr(nbpt) = jg IF (inde.EQ.0) THEN SEGINI,inde ELSE IF (jg.GT.inde(/1)) THEN SEGADJ,inde ENDIF DO j = 1, jg inde(j) = 0 ENDDO ENDIF DO j = 1, nbel1 DO m = 1, nbno1 ib = ipt1.num(m,j) ia = icpr(ib) inde(ia) = j icpr(ib) = ia - 1 ENDDO ENDDO C Fin du travail preparatoire pour le maillage ipt1 ipt2 = imodel.imamod nbno2 = ipt2.num(/1) nbel2 = ipt2.num(/2) c* ipt2 = imodel.imamod = meleme c* nbno2 = ipt2.num(/1) = nbno1 c* nbel2 = ipt2.num(/2) = mlent2.lect(is) C on fabrique le mlenti de correspondance C on dimensionne au nombre d elements de ipt2 = sous-modele is jg = nbel2 SEGINI,mlenti ibon = 0 DO 110 iel2 = 1, nbel2 ia = ipt2.num(1,iel2) ideb = icpr(ia)+1 ifin = icpr(ia+1) IF (ifin.LT.ideb) GOTO 110 DO 111 ib = ideb, ifin iel1 = inde(ib) DO j = 1, nbno1 IF (ipt2.num(j,iel2).NE.ipt1.num(j,iel1)) GOTO 111 ENDDO ibon = ibon + 1 mlenti.lect(iel2) = iel1 GOTO 110 111 CONTINUE 110 CONTINUE IF (ibon .EQ. 0) THEN C Intersection VIDE entre MELEME et IPT1 ixx = 0 ismel(iz,is) = 0 SEGSUP,mlenti ELSE C Intersection NON VIDE entre MELEME et IPT1 IF (ibon.GT.nbel1) THEN C Si on a plus d'elements dans l'intersection que dans ipt1 ! write(ioimp,*) 'REDUAF : Etiquette 11x intersection ?' ENDIF NL3 = NL3 + 1 mlent3.lect(NL3) = mlenti ixx = -1 ismel(iz,is) = NL3 ENDIF C Ajout dans le PRECONDITIONNEMENT : Ajout a la suite IF(mlenti .NE. 0)THEN IPLACE=MOD(NINTSA(ith1),MIN(NTRIPL,max(1,NBESCR)))+1 C PRINT *,'REDUAF_AJOUT',oothrd,IPLACE,meleme,ipt1,mlenti PMAMOD(IPLACE,ith1) = meleme PMAMOH(IPLACE,ith1) = IHO1 PMACHA(IPLACE,ith1) = ipt1 PMACHH(IPLACE,ith1) = IHO2 PMLENT(IPLACE,ith1) = mlenti NINTSA(ith1) = IPLACE ENDIF ENDIF CG write(*,*) ' -',iz,is,ixx,ismel(iz,is) 102 CONTINUE C Sous-zone du mchelm a traiter IF (ixx .NE. 0) THEN DO 105 ia = 1, iz-1 ib = izone(ia,is) IF (ib.EQ.0) GOTO 105 IF (conche(ia)(1:NCONCH).NE.conloc) GOTO 105 DO k = 1, N3 IF (k.NE.4) THEN IF (infche(ia,k).NE.infche(iz,k)) GOTO 105 ENDIF ENDDO izone(iz,is) = ib GOTO 106 105 CONTINUE ISOZM = ISOZM + 1 izone(iz,is) = ISOZM 106 CONTINUE ENDIF CG write(*,*) ' -',iz,is,ixx,izone(iz,is) 101 CONTINUE 100 CONTINUE IF (icpr.NE.0) SEGSUP,icpr IF (inde.NE.0) SEGSUP,inde C --------------------------------- C Construction du MCHELM resultat : C --------------------------------- C Grace au traitement ci-dessus (boucle 105), ISOZM correspond a N1 : N1 = ISOZM L1 = mchelm.titche(/1) N3 = mchelm.infche(/2) CALL oooprl(1) SEGINI,mchel2 mchel2.titche = mchelm.titche mchel2.ifoche = mchelm.ifoche C Pour chaque sous-modele "is", on regroupe les sous-zones du mchelm "iz" C associees (izone(iz,is) > 0) : DO 200 is = 1, NSMOD imodel = kmodel(is) IF (imodel.nefmod.EQ.259) GOTO 200 ipt2 = imodel.imamod nbel2 = mlent2.lect(is) DO 210 iz = 1, NZ in1 = izone(iz,is) IF (in1.LE.0) GOTO 210 mchaml = mchelm.ichaml(iz) n21 = mchaml.ielval(/1) C Cas particulier du mchaml sans composante (on ne fait rien) : IF (n21.EQ.0) GOTO 210 IF (mchel2.imache(in1).EQ.0) THEN CG write(ioimp,*) ' Cas 1 :',mchel2.imache(in1) mchel2.conche(in1) = mchelm.conche(iz) mchel2.imache(in1) = ipt2 DO k = 1, N3 mchel2.infche(in1,k) = mchelm.infche(iz,k) ENDDO n22 = 0 n2 = n22 + n21 SEGINI,mcham2 mchel2.ichaml(in1) = mcham2 ELSE CG write(ioimp,*) ' Cas 2 :',mchel2.imache(in1) mcham2 = mchel2.ichaml(in1) n22 = mcham2.ielval(/1) n2 = n22 + n21 SEGADJ,mcham2 ENDIF *jk148537 ** if (lzsxx(iz,is)) then ** do i = 1,n21 ** mcham2.nomche(n22+i) = mchaml.nomche(i) ** mcham2.typche(n22+i) = mchaml.typche(i) ** mcham2.ielval(n22+i) = mchaml.ielval(i) ** enddo ** goto 210 ** endif mlenti = ismel(iz,is) IF (mlenti.GT.0) mlenti = mlent3.lect(mlenti) CG write(ioimp,*) ' :',iz,is,mlenti,n22,n21,n2 DO i = 1, n21 nomloc = mchaml.nomche(i) iplac = 0 IF (n22.NE.0) THEN ENDIF typloc = mchaml.typche(i) melval = mchaml.ielval(i) nbpi1 = MAX(melval.velche(/1),melval.ielche(/1)) nbel1 = MAX(melval.velche(/2),melval.ielche(/2)) IF (nbel1.GT.1) nbel1 = nbel2 IF (iplac.EQ.0) THEN iplac = n22 + i mcham2.nomche(iplac) = nomloc mcham2.typche(iplac) = typloc IF (typloc.EQ.'REAL*8 ') THEN n1ptel = nbpi1 n1el = nbel2 n2ptel = 0 n2el = 0 ELSE n1ptel = 0 n1el = 0 n2ptel = nbpi1 n2el = nbel2 ENDIF SEGINI,melva2 if (n1ptel.eq.0.and.n2ptel.eq.0) then * write(6,*) 'reduaf melva2 melval ', * > melva2,melval,n1ptel,n1el,n2ptel,n2el endif mcham2.ielval(iplac) = melva2 ELSE C incompatibilite du type de composante entre champs IF (mcham2.typche(iplac).NE.typloc) THEN KERRE = 917 MOTERR(1:4) = nomloc MOTERR(5:21) = typloc MOTERR(22:38) = mcham2.typche(iplac) GOTO 9000 ENDIF melva2 = mcham2.ielval(iplac) * on duplique melva2 au cas ou il soit partage car on va le modifier ** segini,melva3=melva2 ** mcham2.ielval(iplac)=melva3 ENDIF melva2 = mcham2.ielval(iplac) C On ajoute melval a melva2 en tenant compte de l'intersection entre C les maillages (mlenti = 0 si maillage identique, >0 sinon) C "Extension" de melva2 si besoin par rapport a melval (appel a MELEXT) C sera effectuee en prealable de l'addition des valeurs dans MELADD. C si melva2 existait avant, on le duplique avant de le modifier CALL oooho1(melva2,ihmelv) if (ihmelv.ne.ihcour) then segini,melva3=melva2 melva2=melva3 mcham2.ielval(iplac)=melva2 endif IF (KERRE.NE.0) GOTO 9000 ENDDO C 210 CONTINUE 200 CONTINUE C Compactage du champ resultat : C ------------------------------ n1max = n1 n1 = 0 DO 310 i = 1, n1max mcham2 = mchel2.ichaml(i) IF (mcham2.EQ.0) GOTO 310 C on compacte les composantes (s'il y en a bien sur !) n22 = mcham2.ielval(/1) IF (n22.EQ.0) GOTO 312 n2 = 0 DO 311 j = 1, n22 melva2 = mcham2.ielval(j) IF (melva2.EQ.0) GOTO 311 CALL oooho1(melva2,ihmelv) IF(ihmelv .EQ. ihcour)THEN C Reduction seulement pour les SEGMENTS nouveaux ! ENDIF segact melva2 n2 = n2 + 1 mcham2.nomche(n2) = mcham2.nomche(j) mcham2.typche(n2) = mcham2.typche(j) mcham2.ielval(n2) = melva2 311 CONTINUE IF (n2.EQ.0) GOTO 310 IF (n2.NE.n22) SEGADJ,mcham2 312 CONTINUE n1 = n1 + 1 mchel2.conche(n1) = mchel2.conche(i) mchel2.imache(n1) = mchel2.imache(i) mchel2.ichaml(n1) = mcham2 DO j = 1, N3 mchel2.infche(n1,j) = mchel2.infche(i,j) ENDDO 310 CONTINUE IF (n1.NE.n1max) SEGADJ,mchel2 CALL oooprl(0) C Definition du type du MCHAML C typ1 contient le nom du type identifie C ltyp1 la longueur de la chaine de caractere C IF (IERR.NE.0) RETURN C Cas particuliers des modeles de modele (melange) IF(ltyp1.NE.-2 .AND. ltyp1.GT.0 .and. mchel2.titche.eq.' ')THEN IF (ltyp1 .NE. L1 ) THEN L1=ltyp1 SEGADJ, mchel2 ENDIF mchel2.titche=typ1 ENDIF C On sort un champ vide s'il n'y a pas de zone commune : c* IF (n1.EQ.0) THEN c**G if (iimpi.eq.7203) write(ioimp,*) 'N1 = 0 apres traitement' c* KERRE = 21 c* ENDIF 9000 CONTINUE C Destruction des segments de travail devenus inutiles : SEGSUP,izone,ismel,mlent3,mlent2,szsxx 9010 CONTINUE IF (KERRE.NE.0) THEN iret = 0 mchel2 = 0 ENDIF CG if (iimpi.eq.7203) then CG write(ioimp,*) 'Sortie de reduaf',mchel2,kerre CG if (kerre.eq.0) call zpchel(mchel2,1) CG endif C Mise a jour du preconditionnement dans CCPRECO (Nouveau champ mchel2) PRECM2(1,ith1) = mchel2 END
© Cast3M 2003 - Tous droits réservés.
Mentions légales