Télécharger prom.eso

Retour à la liste

Numérotation des lignes :

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

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