Télécharger prom.eso

Retour à la liste

Numérotation des lignes :

  1. C PROM SOURCE FD218221 14/01/31 21:15:12 7936
  2.  
  3. SUBROUTINE PROM(IPMODE,IPCHE,IPCHT,ISUP,IPOUT)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. EXTERNAL SHAPE
  7. *-------------------------------------------------------------------- *
  8. * *
  9. * Sous-programme appelé par l'opérateur PROI *
  10. * Projection d'un chamelem sur le support d un modele par *
  11. * minimisation de l integrale (Ue-Us)**2 ds sur chaque element *
  12. * Ue evalue a partir des fonctions de formes du modele associe *
  13. * au champ d origine si il est fourni ( sinon a partir des *
  14. * fonctions de transformation geometriques de l element) *
  15. * *
  16. * L integrale est exprime à partir des fonctions de forme associees*
  17. * au modèle support recepteur *
  18. *---------------------------------------------------------------------*
  19. * *
  20. * ENTREE : *
  21. * IPMODE modele sur lequel on va projeter *
  22. * IPCHE MCHAML de caracteristique (si le maodele comprends *
  23. * des coques) *
  24. * IPCHT MCHAML a projeter (obligatoirement defini aux noeuds)
  25. * ISUP remis a 1 pour le moment (on travaille aux noeuds)*
  26. * *
  27. * SORTIE : *
  28. * IPOUT MCHAML resultat *
  29. * *
  30. *---------------------------------------------------------------------*
  31. *
  32. -INC CCOPTIO
  33. -INC SMCOORD
  34. -INC SMMODEL
  35. -INC SMCHAML
  36. -INC SMELEME
  37. -INC SMCHPOI
  38. -INC SMLCHPO
  39. -INC SMLMOTS
  40. -INC TMTRAV
  41. -INC SMINTE
  42. -INC CCREEL
  43. *
  44. PARAMETER ( O1=1.D0,O2=2.D0)
  45. character*4 nomm
  46. segment snomip
  47. character*4 nomip(0)
  48. endsegment
  49. * SEGMENT MDATA
  50. * INTEGER IPCA(*),IPEXI(*)
  51. * ENDSEGMENT
  52. segment ipcb(nsschp)
  53. segment icpr(2,xcoor(/1)/idim+1)
  54. segment icpr1(xcoor(/1)/idim+1)
  55. C
  56. segment info
  57. integer infell(JG)
  58. endsegment
  59. C
  60. segment sicoq
  61. integer icocoq(nbsous),imfr(nbsous),iiele(nbsous)
  62. integer iminte(nbsous),imint1(nbsous),ibnap(nbsous)
  63. integer inbg(nbsous),ibbas(nbsous)
  64. endsegment
  65. C
  66. segment WRK5
  67. REAL*8 W(200),A(NBN1,NBN1),B(NBN1,NBN1),XPU(3),RHS(NBN1)
  68. REAL*8 SHPXXX(6,NBNN,200),XX(3,200),VAL(NBN1),POIDS(200)
  69. INTEGER NDD
  70. endsegment
  71. C
  72. SEGMENT WRK4
  73. REAL*8 BPSS(3,3),XEL(3,NBN1),XU(3),SHP(6,NBN1),S(6)
  74. REAL*8 TXR(3,3,NBN1),XE(3,NBG)
  75. ENDSEGMENT
  76. SEGMENT VECT
  77. REAL*8 VEC1(IDIM)
  78. REAL*8 VEC2(IDIM)
  79. REAL*8 VECN(IDIM,ibbas(isous))
  80. ENDSEGMENT
  81.  
  82. C= Logique a TRUE si mode de calcul est axisymetrique 1D/2D ou Fourier
  83. LOGICAL LOGAXI
  84. C= Quelques constantes (2.Pi et 4.Pi)
  85. PARAMETER (X2Pi =6.283185307179586476925286766559D0)
  86. PARAMETER (X4Pi=12.566370614359172953850573533118D0)
  87.  
  88. C on remet isup a 1 pour le moment
  89. ISUP=1
  90. MMODEL=IPMODE
  91. SEGACT,MMODEL
  92. NSOUS=KMODEL(/1)
  93. C
  94. ndec=0
  95. call lirent(ndec,0,iret)
  96. if(ndec.eq.0) ndec = 3
  97. *
  98. NBNOE = XCOOR(/1)/(IDIM+1)
  99.  
  100. if(ipche.NE.0) then
  101. MCHELM=IPCHE
  102. segact MCHELM
  103. NCHAM=ICHAML(/1)
  104. endif
  105. *
  106. * Création du maillage provisoire support pour interpolation
  107. *
  108. NBSOUS = NSOUS
  109. NBREF = 0
  110. NBELEM = 0
  111. NBNN = 0
  112. C ipt2 contiendra les super couche ou points de gauss
  113. SEGINI IPT2
  114. SEGINI SICOQ
  115. c listmots des phases
  116. ilphmo = -1
  117. jgn = 8
  118. jgm = nsous
  119. segini mlmots
  120. ilphmo = mlmots
  121. jgm = 1
  122. *
  123. * Boucle sur l'ensemble des sous zones du modeles sur lequel on va
  124. * projetter pour initialisations
  125. nponou=0
  126. npcoq =0
  127. DO 30 ISOUS = 1, NSOUS
  128. *
  129. IMODEL = KMODEL(ISOUS)
  130. SEGACT,IMODEL
  131. MELEME = IMAMOD
  132. SEGACT MELEME
  133. MELE = NEFMOD
  134. NBELE1= NUM(/2)
  135. NBN1 = NUM(/1)
  136. ibnap(isous)=1
  137. C
  138. C COQUE INTEGREE OU PAS ?
  139. C
  140. IF(INFMOD(/1).NE.0)THEN
  141. NPINT=INFMOD(1)
  142. if( npint.ne.0)IBNAP(isous) =NPINT
  143. ELSE
  144. NPINT=0
  145. C ibnap(isous)=1
  146. ENDIF
  147.  
  148. if( infmod(/1).lt.2+isup) then
  149. call elquoi(mele,npint,isup,ipinf,imodel)
  150. if (ierr.ne.0) then
  151. segdes imodel,mmodel
  152. segsup mchelm
  153. iret=0
  154. return
  155. endif
  156. info=ipinf
  157. mfr =infell(13)
  158. iele =infell(14)
  159. minte=infell(11)
  160. minte1=infell(12)
  161. nbg=infell(6)
  162. segsup info
  163. else
  164. minte=infmod(isup+2)
  165. minte1=infMOD(12)
  166. nbg=infele(6)
  167. mfr =infele(13)
  168. iele =infele(14)
  169. endif
  170. imfr(isous)=mfr
  171. iiele(isous)=iele
  172. iminte(isous)=minte
  173. imint1(isous)=minte1
  174. inbg(isous)=nbg
  175. * write(6,*) 'mfr iele mele minte minte1 nbg',mfr,iele,
  176. * & mele,minte,minte1,nbg
  177. NBSOUS = 0
  178. NBREF = 0
  179. NBELEM = NBELE1
  180.  
  181. IF(mfr.eq.3.or.mfr.eq.5.or.mfr.eq.9) THEN
  182. C------------------coques
  183. IF(IPCHE.eq.0) THEN
  184. C pour le coques il faut les caracteristiques
  185. call erreur(404)
  186. segdes meleme ,imodel,mmodel
  187. * segsup info
  188. return
  189. ENDIF
  190. C
  191. icocoq(isous)=1
  192. nbnn = ndec*ndec
  193. if(iele.lt.4) nbnn=ndec
  194. ibbas(isous)= nbnn
  195. if(idim.eq.3.or.(iele.lt.4.and.idim.eq.2)) then
  196. ibnap(isous)=ndec
  197. nbnn= nbnn*ndec
  198. inbg(isous)=nbnn
  199. endif
  200.  
  201. ELSE
  202. C--------------- massifs
  203.  
  204. IF (mfr.EQ.1) THEN
  205. IF (IDIM.EQ.3) THEN
  206. NBNN=ndec*ndec*ndec
  207. ELSE IF (IDIM.EQ.2) THEN
  208. NBNN=ndec*ndec
  209. C* ELSE IF (IDIM.EQ.1) THEN
  210. ELSE
  211. NBNN=ndec
  212. ENDIF
  213. ELSE
  214. NBNN=ndec*ndec
  215. ENDIF
  216.  
  217. C petite rectif utile pour les elements speciaux
  218. if(iele.eq.32) then
  219. C les polygones
  220. nbnn =nbnn*nbn1
  221. elseif(iele.eq.25.or.iele.eq.26) then
  222. C les pyramides
  223. nbnn = nbnn*4
  224. endif
  225. C
  226. inbg(isous)= NBNN
  227. ENDIF
  228. *
  229. SEGINI IPT1
  230. IPT1.ITYPEL = 28
  231. ipt2.lisous(isous) = ipt1
  232. segdes ipt1
  233. nponou=nponou+nbnn*nbelem
  234.  
  235. if (isous.eq.1) then
  236. mots(1) = conmod(17:24)
  237. else
  238. do ipl = 1,jgm
  239. if (mots(ipl).eq.conmod(17:24)) goto 27
  240. enddo
  241. jgm = jgm + 1
  242. mots(jgm) = conmod(17:24)
  243. 27 continue
  244. endif
  245.  
  246. segdes meleme,imodel
  247. * segsup info
  248. 30 CONTINUE
  249. segadj mlmots
  250. C-----------------------------------------------------------
  251. C
  252. C
  253. C on va fabriquer des points support provisoires a l interieur
  254. C de chaque element des maillages du modele recepteur
  255. C
  256. C
  257. C----------------------------------------------------------
  258. C write(6,*) ' nombre de points crees ' ,nponou
  259. C (DIP) ajustement de MCOORD pour les noeuds provisoires
  260. C il sera remis a sa taille initiale a la fin du programme
  261. NBPTS = NBNOE + nponou
  262. SEGADJ MCOORD
  263. *-----------------------------------------------------------
  264. * Boucle sur l'ensemble des sous zones du modeles
  265. *------------------------------------------------------------
  266. inupo =NBNOE
  267. DO 100 ISOUS = 1,NSOUS
  268. IMODEL = KMODEL(ISOUS)
  269. segact imodel
  270. MELEME = IMAMOD
  271. segact meleme
  272. nbele1= num(/2)
  273. nbn1 = num(/1)
  274. ipt1=ipt2.lisous(isous)
  275. segact ipt1*mod
  276. NNB = ibbas(isous)
  277. NBG = inbg(isous)
  278. SEGINI WRK4
  279. MELVA1 = 0
  280. MELVA2 = 0
  281. IF(icocoq(isous).ne.0) THEN
  282. C -------------------------------la sz est une coque ------
  283.  
  284. NUCHA = 0
  285. DO 15, NUCH = 1, NCHAM
  286. *
  287. IF ( CONCHE(NUCH).EQ.CONMOD.AND.
  288. & IMACHE(NUCH).EQ.IMAMOD) NUCHA = NUCH
  289. *
  290. 15 CONTINUE
  291. *
  292. IF (NUCHA.NE.0) THEN
  293. MCHAML=ICHAML(NUCHA)
  294. SEGACT,MCHAML
  295. *
  296. MELVA1 = 0
  297. MELVA2 = 0
  298. NCOMP = IELVAL(/1)
  299. DO 20, I = 1, NCOMP
  300. IF (NOMCHE(I)(1:4).EQ.'EPAI') THEN
  301. MELVA1 = IELVAL(I)
  302. SEGACT, MELVA1
  303. ELSEIF (NOMCHE(I)(1:4).EQ.'EXCE') THEN
  304. MELVA2 = IELVAL(I)
  305. SEGACT, MELVA2
  306. ENDIF
  307. 20 CONTINUE
  308. ENDIF
  309. ENDIF
  310. C
  311. MELE = NEFMOD
  312. c coque integree ou pas ?
  313. if(infmod(/1).ne.0)then
  314. npint=infmod(1)
  315. else
  316. npint=0
  317. endif
  318. C
  319. nbnap=ibnap(isous)
  320. mfr =imfr(isous)
  321. iele =iiele(isous)
  322. minte=iminte(isous)
  323. minte1=imint1(isous)
  324. nbg = inbg(isous)
  325. C
  326. segact minte
  327. if(mfr.eq.5) segact minte1
  328.  
  329. nnb = ibbas(isous)
  330. C
  331. ipo=0
  332. *-------------------------------------------------------------------
  333. * Pavage de l element de reference
  334. *-------------------------------------------------------------------
  335. NBNN = inbg(isous)
  336. segini WRK5
  337. iwkr= wrk5
  338. iele=iiele(isous)
  339. CALL PAVPOI(WRK5,iele,ndec,nbn1)
  340. segini vect
  341. *-------------------------------------------------------------------
  342. * Boucle sur les elements de la sous zone du modele
  343. *-------------------------------------------------------------------
  344.  
  345. DO IEL=1,NBELE1
  346. CALL DOXE(XCOOR,IDIM,NBN1,NUM,IEL,XEL)
  347.  
  348. C si c est une coque pouvant etre excentree
  349. C recherche des normales aux noeuds
  350. if(icocoq(isous).ne.0) then
  351. C
  352. IF(IELE.EQ.4) THEN
  353. C --------coq3 dkt dst ----------------
  354. do 98 i=1,idim
  355. VEC1(i)=xel(i,2)-xel(i,1)
  356. 98 VEC2(i)=xel(i,3)-xel(i,1)
  357. vecn(1,1) = vec1(2)*vec2(3)-vec1(3)*vec2(2)
  358. vecn(2,1) = vec1(3)*vec2(1)-vec1(1)*vec2(3)
  359. vecn(3,1) = vec1(1)*vec2(2)-vec1(2)*vec2(1)
  360. vnor =sqrt(vecn(1,1)*vecn(1,1)+vecn(2,1)*vecn(2,1)
  361. & +vecn(3,1)*vecn(3,1))
  362. do 91 i=1,idim
  363. 91 vecn(i,1)=vecn(i,1)/vnor
  364. do 99 i=1,idim
  365. do 99 j=2,ibbas(isous)
  366. 99 vecn(i,j)=vecn(i,1)
  367. ELSEIF (IELE.EQ.8) THEN
  368. C --------------- coq4 ----------------
  369. call cq4loc(xel,xe,BPSS,irrt,1)
  370. C ici c est une estimation du plan moyen
  371. * write(6,fmt='(3E12.5)') ((bpss(i,j),i=1,3),j=1,3)
  372. do i=1,3
  373. do ip=1,ibbas(isous)
  374. vecn(i,ip)=bpss(3,i)
  375. enddo
  376. enddo
  377. ELSEIF(IELE.eq.2) THEN
  378. C ---------------- coq2 --------------
  379. vnor=0.D0
  380. do 92 i=1,idim
  381. VEC1(i)=xel(i,2)-xel(i,1)
  382. 92 vnor=vnor + vec1(i)*vec1(i)
  383. vnor = sqrt(vnor)
  384. vecn(1,1)=-vec1(2)/vnor
  385. vecn(2,1)=vec1(1)/vnor
  386. vecn(1,2)=vecn(1,1)
  387. vecn(2,2)=vecn(2,1)
  388. ELSEIF (iele.eq.5.or.iele.eq.10) THEN
  389. C preparation pour coques integrees
  390. SURFX=XZERO
  391. SURFY=XZERO
  392. SURFZ=XZERO
  393. ENDIF
  394. endif
  395. C ------------------- boucle sur les points d integration
  396. inuo=inupo
  397. C write(6,*) ' ndd et inupo ' ,ndd,inupo
  398. do igau=1,ndd
  399. C
  400. do j=1,idim
  401. xpu(j)=xx(j,igau)
  402. enddo
  403. C write(6,*)(xpu(J),j=1,idim)
  404. C
  405. if(mele.GE.111.and.mele.le.122) then
  406. C les plolygones
  407. call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  408. else
  409. CALL SHAPE(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  410. endif
  411.  
  412. C PROJECTION DANS L ELEMENT REEL
  413. DO 501 IK = 1,IDIM
  414. XU(IK)= XZERO
  415. DO 501 JK = 1,nbn1
  416. XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK))
  417. 501 CONTINUE
  418. C
  419. if ( mfr.eq.5) then
  420. C coques integrees ds l epaisseur
  421. DO I=1,6
  422. S(I)=XZERO
  423. CONTINUE
  424. ENDDO
  425. DO INOE=1,NBN1
  426. S(1)=S(1)+SHP(2,INOE)*XEL(2,INOE)
  427. S(2)=S(2)+SHP(3,INOE)*XEL(3,INOE)
  428. S(3)=S(3)+SHP(3,INOE)*XEL(2,INOE)
  429. S(4)=S(4)+SHP(2,INOE)*XEL(3,INOE)
  430. S(5)=S(5)+SHP(3,INOE)*XEL(1,INOE)
  431. S(6)=S(6)+SHP(2,INOE)*XEL(1,INOE)
  432. ENDDO
  433. SURFX=S(1)*S(2)-S(3)*S(4)
  434. SURFY=S(4)*S(5)-S(2)*S(6)
  435. SURFZ=S(6)*S(3)-S(5)*S(1)
  436. SURF=SQRT(SURFX**2+SURFY**2+SURFZ**2)
  437. C write(6,*) ' normales au pt igau du plan de base '
  438. vecn(1,igau)=surfx/surf
  439. vecn(2,igau)=surfy/surf
  440. vecn(3,igau)=surfz/surf
  441. endif
  442. C
  443. IF(ICOCOQ(ISOUS).EQ.0) THEN
  444. C elements standards
  445. inupo=inupo+1
  446. ipt1.num(igau,iel)=inupo
  447. do i=1,idim
  448. xcoor((inupo-1)*(idim+1)+i)=XU(i)
  449. enddo
  450. C write(6,3335) (xu(i),i=1,idim)
  451. C3335 format(3e12.5)
  452. ELSE
  453. C traitement special des coques
  454. if(melva1.ne.0) then
  455. ibmn =MIN(iel,melva1.velche(/2))
  456. igmn =MIN(igau,melva1.velche(/1))
  457. epai = melva1.velche(igmn,ibmn)
  458. else
  459. epai =xzero
  460. endif
  461. C
  462. if(melva2.ne.0) then
  463. ibmn =MIN(iel,melva2.velche(/2))
  464. igmn =MIN(igau,melva2.velche(/1))
  465. exce = melva2.velche(igmn,ibmn)
  466. else
  467. exce =xzero
  468. endif
  469.  
  470. dp = epai/ibnap(isous)
  471. dnor = exce -(epai + dp)/2.D0
  472. C
  473. do inap=1,ibnap(isous)
  474. dnor=dnor+dp
  475. ipos = (inap-1)*ibbas(isous)+igau
  476. ipt1.num(ipos,iel)=inuo+ipos
  477. do i=1,idim
  478. xcoor((inuo+ipos-1)*(idim+1)+i)= XU(i)+dnor*vecn(i,igau)
  479. enddo
  480. enddo
  481. inupo=inupo+ibnap(isous)
  482. ENDIF
  483.  
  484. C ---------------fin de la boucle sur les points de d integration
  485. enddo
  486.  
  487. C fin de la boucle sur les elements
  488. enddo
  489. segsup wrk5,vect
  490. segdes minte
  491. if(mfr.eq.5) segdes minte1
  492. if(icocoq(isous).ne.0) SEGDES,MCHAML
  493. if (melva1.ne.0) segdes melva1
  494. if (melva2.ne.0) segdes melva2
  495. segdes meleme,imodel
  496. segdes ipt1
  497. SEGSUP WRK4
  498. * segsup info
  499. *
  500. 100 CONTINUE
  501. C fin de creation des noeuds supports provisoires
  502. 2003 format(3I4,2X,4e12.5,I4)
  503. C--------------//////////////////////////////////------------------
  504. C maintenant in faut un meleme poi1 de tous les points sur lesquels
  505. C on doit interpoler une valeur
  506. segini icpr1
  507. ipt5=ipt2
  508. nbelem=0
  509. do io=1,max(1,ipt2.lisous(/1))
  510. if (ipt2.lisous(/1).ne.0) ipt5=ipt2.lisous(io)
  511. segact ipt5
  512. do ip=1,ipt5.num(/1)
  513. do il=1,ipt5.num(/2)
  514. ipt=ipt5.num(ip,il)
  515. if (icpr1(ipt).eq.0) then
  516. nbelem=nbelem+1
  517. icpr1(ipt)=nbelem
  518. endif
  519. enddo
  520. enddo
  521. ** segdes ipt5
  522. enddo
  523. ** segdes ipt2
  524. nbnn=1
  525. nbsous=0
  526. nbref=0
  527. ** write (6,*) ' nbelem ',nbelem
  528. segini ipt4
  529. ipgeom=ipt4
  530. nbelem=0
  531. do i=1,icpr1(/1)
  532. if (icpr1(i).ne.0) then
  533. nbelem=nbelem+1
  534. ** write (6,*) ' i nbelem icpr ',i,nbelem,icpr1(i)
  535. ipt4.num(1,nbelem)=i
  536. endif
  537. enddo
  538. segsup icpr1
  539. C----------------------------------------------------------------
  540. C----------------------------------------------------------------
  541. *-------------on est pret a faire l interpolation sur ipgeom ----
  542. C
  543. *write(6,*) ' maillage apres changer poi1 ' ,ipgeom
  544. C call ecmail(ipgeom)
  545. isort=0
  546. CALL PRO2(IPGEOM,IPCHT,isort,IPOUT,ilphmo)
  547. if (ierr.ne.0) return
  548. C----------------------------------------------------------------
  549. C----------------------------------------------------------------
  550. C
  551. C Definition de "constantes" selon le mode de calcul
  552. C Modes de calcul 1D/2D axisymetrique
  553. IF (IFOMOD.EQ.0.OR.IFOMOD.EQ.4) THEN
  554. LOGAXI=.TRUE.
  555. CteAxi=X2Pi
  556. C Mode de calcul 2D Fourier
  557. ELSE IF (IFOMOD.EQ.1) THEN
  558. LOGAXI=.TRUE.
  559. IF (NIFOUR.EQ.0) THEN
  560. CteAxi=X2Pi
  561. ELSE
  562. CteAxi=XPi
  563. ENDIF
  564. ELSE
  565. LOGAXI=.FALSE.
  566. ENDIF
  567. C----------------------------------------------------------------
  568. C --- il faut maintenant reconstituer les MCHAML
  569. C --- a partir du chpo construit sur ipgeom ----------------
  570. mlchpo = ipout
  571. segact mlchpo
  572. mchel5= ipcht
  573. segact mchel5
  574. N1=NSOUS
  575. L1=12
  576. N3=6
  577. SEGINI MCHEL1
  578. * MCHEL1.TITCHE='SCALAIRE'
  579. MCHEL1.TITCHE = mchel5.titche
  580. MCHEL1.IFOCHE=IFOUR
  581. segdes mchel5
  582. *//////////////////////////////
  583.  
  584. * boucle phases modele
  585. DO 7000 IPHAS = 1,JGM
  586. * write(6,*) 'iphas ',iphas
  587. segini icpr,snomip
  588. * mdata=isort
  589. * segact mdata
  590. mchpoi=ichpoi(iphas)
  591. segact mchpoi
  592. nsschp=ipchp(/1)
  593. segini ipcb
  594.  
  595. C ip=0
  596. do i=1,nsschp
  597. * mchpoi=ipca(i)
  598. C call ecchpo(mchpoi)
  599. * segact mchpoi
  600. msoupo=ipchp(i)
  601. ipcb(i)=msoupo
  602. segact msoupo
  603. mpoval=ipoval
  604. segact mpoval
  605. nc = nocomp(/2)
  606. C
  607. C tableau des general des composantes
  608. do ic = 1,nc
  609. nomm= nocomp(ic)
  610. if(nomip(/2).eq.0) then
  611. nomip(**)=nomm
  612. else
  613. do k=1,nomip(/2)
  614. if(nomm.eq.nomip(k)) goto 4322
  615. enddo
  616. nomip(**)=nomm
  617. 4322 continue
  618. endif
  619. enddo
  620. * write(6,*) (nomip(l),l=1,nomip(/2))
  621. * reperage de la position ses points dans le chpo
  622. ipt5=igeoc
  623. segact ipt5
  624. * ip=0
  625. do j=1,ipt5.num(/2)
  626. * ip=ip+1
  627. icpr(1,ipt5.num(1,j))=j
  628. icpr(2,ipt5.num(1,j))= i
  629. enddo
  630. enddo
  631. C i=icpr(2,k) point k venant du msoupo du msoupo i dans ipcb
  632.  
  633. inupo=0
  634. segact ipt2
  635. DO 2100 isous = 1,nsous
  636. imodel= KMODEL(isous)
  637. segact imodel
  638. if (conmod(17:24).ne.mots(iphas)) goto 2100
  639. meleme = imamod
  640. segact meleme
  641. idi2=idim
  642. ipt5=ipt2.lisous(isous)
  643. segact ipt5
  644. if(icocoq(isous).ne.0) idi2 = idim-1
  645. mfr=imfr(isous)
  646. * création du nouveau segment MCHAML
  647. N2=nomip(/2)
  648. SEGINI MCHAML
  649. MCHEL1.IMACHE(isous)=MELEME
  650. MCHEL1.ICHAML(isous)=MCHAML
  651. MCHEL1.CONCHE(isous)=CONMOD
  652. MCHEL1.INFCHE(isous,1)=0
  653. MCHEL1.INFCHE(isous,2)=0
  654. MCHEL1.INFCHE(isous,3)=NIFOUR
  655. MCHEL1.INFCHE(isous,4)=0
  656. MCHEL1.INFCHE(isous,5)=0
  657. MCHEL1.INFCHE(isous,6)=1
  658. *
  659. N1EL=NUM(/2)
  660. N2PTEL=0
  661. N2EL=0
  662. C--------------------------------------------------------
  663. C boucle sur les composantes
  664. C--------------------------------------------------------
  665. NBNN=num(/1)
  666. N1PTEL= NBNN
  667. NBN1=NBNN
  668.  
  669. C
  670. SEGINI WRK5,wrk4
  671. iele = iiele(isous)
  672. iwrk= wrk5
  673. CALL PAVPOI(WRK5,iele,ndec,nbnn)
  674. C ----- remplissage de shpxxx ( shptot)
  675. DO IGAU=1,NDD
  676. do j=1,idim
  677. xpu(j)=xx(j,igau)
  678. enddo
  679. if(mele.GE.111.and.mele.le.122) then
  680. call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  681. else
  682. call shape(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  683. endif
  684. do id=1,6
  685. do no=1,nbnn
  686. shpxxx(id,no,igau)=shp(id,no)
  687. enddo
  688. enddo
  689. ENDDO
  690. N1PTEL=NBNN
  691. C ----------- creation des melvals recepteurs -----
  692. DO icomp=1,NOMIP(/2)
  693. SEGINI MELVAL
  694. NOMCHE(icomp)=nomip(icomp)
  695. TYPCHE(icomp)='REAL*8'
  696. IELVAL(icomp)=MELVAL
  697. ENDDO
  698. C--------------------------------------------------
  699. C-------- boucle sur les elements ----------
  700. AXIS=O1
  701. DO 2210 iel=1,num(/2)
  702. CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XEL)
  703. CALL ZERO(A,NBNN,NBNN)
  704. C SSOM1=XZERO
  705. C
  706. C
  707. DO inap=1,ibnap(isous)
  708. DO 40 IGAU=1,ndd
  709. do i=1,idim
  710. xpu(i)=xx(i,igau)
  711. enddo
  712. if(mele.GE.111.and.mele.le.122) then
  713. C les plolygones
  714. call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  715. else
  716. call shape(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  717. endif
  718. CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC)
  719. C Modes de calcul 1D/2D axisymetrique et 2D Fourier
  720. IF (LOGAXI) THEN
  721. CALL DISTRR(XEL,SHP,NBNN,RR)
  722. AXIS=CteAxi*RR
  723. C Mode de calcul 1D spherique
  724. ELSE IF (IFOMOD.EQ.5) THEN
  725. CALL DISTRR(XEL,SHP,NBNN,RR)
  726. AXIS=X4Pi*RR*RR
  727. ENDIF
  728. C
  729. C on veut calculer la fonction au point d integration projete
  730. C PROJECTION
  731. DO IK = 1,IDIM
  732. XU(IK)= XZERO
  733. DO JK = 1,NBNN
  734. XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK))
  735. ENDDO
  736. ENDDO
  737. C
  738. CPZA=poids(igau)*DJAC*AXIS
  739. DO 50 INO1=1,NBNN
  740. DO 60 INO2=1,NBNN
  741. Q = SHPXXX(1,INO1,IGAU)*SHPXXX(1,INO2,IGAU)*CPZA
  742. CCC Q = SHPXXX(1,INO1,IGAU)*SHPXXX(1,INO2,IGAU)
  743. CCC 1 *poids(igau)*DJAC*AXIS
  744. A(INO1,INO2)=A(INO1,INO2)+Q
  745. 60 CONTINUE
  746. 50 CONTINUE
  747. 40 CONTINUE
  748. enddo
  749. C
  750. C ---- fin preparation matrice 'masse ' sur element courant
  751. C do i=1,NBNN
  752. C do j=1,NBNN garder eventuellement pour calculer
  753. C B(i,j)=A(i,j) la difference sur les integrales
  754. C enddo
  755. C enddo
  756. C
  757. call invalm(A,NBNN,NBNN,KERRE,1.D-10)
  758. C ------------------ boucle sur les composantes
  759. DO 2220 icomp=1,nomip(/2)
  760. idebco = inupo
  761. do i=1,nbnn
  762. RHS(i) =XZERO
  763. enddo
  764. C --------- boucle sur les points d integration
  765. C
  766. AXIS = O1
  767. DO 42 inap=1,ibnap(isous)
  768. DO 41 IGAU=1,ndd
  769. do i=1,idim
  770. xpu(i)=xx(i,igau)
  771. enddo
  772. CALL SHAPE(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  773. CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC)
  774. C Modes de calcul 1D/2D axisymetrique et 2D Fourier
  775. IF (LOGAXI) THEN
  776. CALL DISTRR(XEL,SHP,NBNN,RR)
  777. AXIS=CteAxi*RR
  778. C Mode de calcul 1D spherique
  779. ELSE IF (IFOMOD.EQ.5) THEN
  780. CALL DISTRR(XEL,SHP,NBNN,RR)
  781. AXIS=X4Pi*RR*RR
  782. ENDIF
  783.  
  784. ipop=(inap-1)*ndd+igau
  785. inupo=ipt5.num(ipop,iel)
  786.  
  787. jh=icpr (1,inupo)
  788. ipa=icpr(2,inupo)
  789. C il faut verifier si vpocha existe pour ce point
  790. msoupo = ipcb(ipa)
  791. mpoval=ipoval
  792. do l=1,nocomp(/2)
  793. if(nocomp(l).eq.nomip(icomp)) then
  794. vvv=vpocha(jh,icomp)
  795. goto 4555
  796. endif
  797. enddo
  798. vvv=XZERO
  799. 4555 continue
  800. C write(6,7777) inupo,jh,poids(igau),vvv
  801. C
  802. melval=ielval(icomp)
  803. CPZA=vvv*poids(igau)*DJAC*AXIS
  804. DO INO1 = 1,NBNN
  805. RHS(ino1)=RHS(ino1)+CPZA*SHPXXX(1,ino1,igau)
  806. CCC RHS(ino1)=RHS(ino1)+vvv*poids(igau)*
  807. CCC 1 DJAC*SHPXXX(1,INO1,IGAU)*AXIS
  808. ENDDO
  809. 41 continue
  810. 42 continue
  811.  
  812. 7777 format(2I4,2X,2E12.5)
  813. C SSOM1 =SSOM1+SURF*poids(igau)*vvv
  814. C
  815. do i=1,NBNN
  816. val(i)=XZERO
  817. do j=1,NBNN
  818. val(i)=val(i)+A(i,j)*RHS(J)
  819. enddo
  820. enddo
  821. C
  822. do I=1,NBNN
  823. velche(i,iel)=val(i)
  824. enddo
  825. C
  826. * if(icomp.lt.nomip(/2)) inupo =idebco
  827. C ---- fin de la boucle sur les composantes
  828. 2220 continue
  829. C fin de la boucle sur les elements
  830. 2210 continue
  831. C
  832. do i=1,ielval(/1)
  833. ipc1 = ielval(i)
  834. call comred(ipc1)
  835. melval = ipc1
  836. segdes melval
  837. enddo
  838. segdes meleme,mchaml
  839.  
  840. C---------------- fin du traitement des massifs
  841. segsup wrk5,wrk4,ipt5
  842.  
  843. C segsup ipt3
  844. segdes imodel,meleme
  845. C fin de la boucle sur les sous zones du modele
  846. 2100 continue
  847. C destruction des chpo intermediaires
  848. do i=1,ipcb(/1)
  849. msoupo=ipcb(i)
  850. cgf ipt5= igeoc (correction 7284)
  851. mpoval= ipoval
  852. * mchpoi=ipca(i)
  853. segsup mpoval,msoupo
  854. cgf segsup ipt5 (correction 7284)
  855. enddo
  856. segsup mchpoi
  857. segsup ipcb,snomip
  858. segsup icpr
  859. 7000 CONTINUE
  860. C (fdp) re-ajustement de MCOORD a sa taille initiale
  861. NBPTS = NBNOE
  862. SEGADJ MCOORD
  863. C retrait des maillages temporaires du pre-conditionnement
  864. c (leurs numero de noeuds depasse la taille de MCOOR)
  865. call crech1b
  866. segdes mchel1,mmodel
  867. C (fdp) suppression du maillage temporaire IPT1
  868. ipt1=ipgeom
  869. segsup ipt1
  870. segact ipt2
  871. do io=1,ipt2.lisous(/1)
  872. ipt5=ipt2.lisous(io)
  873. * segsup ipt5
  874. enddo
  875. segsup ipt2
  876. segsup ipt4
  877. if(ipche.ne.0) then
  878. mchelm=ipche
  879. segdes mchelm
  880. endif
  881. IPOUT=MCHEL1
  882. segsup sicoq
  883.  
  884. return
  885. END
  886.  
  887.  
  888.  
  889.  
  890.  
  891.  

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