36
37
38
40 USE my_alloc_mod
41
42
43
44#include "implicit_f.inc"
45#include "comlock.inc"
46
47
48
49#include "com01_c.inc"
50#include "com04_c.inc"
51#include "param_c.inc"
52#include "remesh_c.inc"
53#include "scr17_c.inc"
54#include "scr18_c.inc"
55
56
57
58 INTEGER IXC(NIXC,*), IPARTC(*), IXTG(NIXTG,*), IPARTTG(*),
59 . IPART(LIPART1,*), SH4TREE(KSH4TREE,*), SH3TREE(KSH3TREE,*)
60 INTEGER ,INTENT(IN) :: NODADT_THERM
61 INTEGER ,INTENT(IN) :: ITHERM_FE
62 my_real a(3,*), stifn(*), ar(3,*), stifr(*), x(3,*),
63 . stcont(*), fthe(*),condn(*)
64
65
66
67 INTEGER KN, KN1, KN2, KN3, KN4
68 INTEGER N, NN, LEVEL, IP, NLEV
69 INTEGER SON,M1,M2,M3,M4,MC,N1,N2,N3,N4,J,K
70 INTEGER I,LLNOD,
71 . LE,LELT,LEV,NE,LELT1,LELT2,
72 . NI,LL
73 INTEGER, DIMENSION(:), ALLOCATABLE :: LNOD
74 INTEGER, DIMENSION(:), ALLOCATABLE :: NELT
75 INTEGER, DIMENSION(:), ALLOCATABLE :: LKINNOD
77 . a1,a2,a3,a4,ac,
78 . phi,facm,faci,r,s
79 my_real,
DIMENSION(:),
ALLOCATABLE :: rnod
80 my_real,
DIMENSION(:),
ALLOCATABLE :: snod
81
82 CALL my_alloc(lnod,numnod)
83 CALL my_alloc(nelt,2*(4**levelmax))
84 CALL my_alloc(lkinnod,numnod)
85 CALL my_alloc(rnod,numnod)
86 CALL my_alloc(snod,numnod)
87
88 lkinnod=0
89 DO level=levelmax-1,0,-1
90
93
94 son=sh4tree(2,n)
95
96 n1=ixc(2,n)
97 n2=ixc(3,n)
98 n3=ixc(4,n)
99 n4=ixc(5,n)
100
101 mc=ixc(4,son)
102 DO j=1,3
103 ac= fourth*a(j,mc)
104 a(j,n1)=a(j,n1)+ac
105 a(j,n2)=a(j,n2)+ac
106 a(j,n3)=a(j,n3)+ac
107 a(j,n4)=a(j,n4)+ac
108 END DO
109 ac=fourth*stifn(mc)
110 stifn(n1)=stifn(n1)+ac
111 stifn(n2)=stifn(n2)+ac
112 stifn(n3)=stifn(n3)+ac
113 stifn(n4)=stifn(n4)+ac
114 IF(istatcnd/=0)THEN
115 ac=fourth*stcont(mc)
116 stcont(n1)=stcont(n1)+ac
117 stcont(n2)=stcont(n2)+ac
118 stcont(n3)=stcont(n3)+ac
119 stcont(n4)=stcont(n4)+ac
120 END IF
121
122 DO j=1,3
123 ac= fourth*ar(j,mc)
124 ar(j,n1)=ar(j,n1)+ac
125 ar(j,n2)=ar(j,n2)+ac
126 ar(j,n3)=ar(j,n3)+ac
127 ar(j,n4)=ar(j,n4)+ac
128 END DO
129 ac=fourth*stifr(mc)
130 stifr(n1)=stifr(n1)+ac
131 stifr(n2)=stifr(n2)+ac
132 stifr(n3)=stifr(n3)+ac
133 stifr(n4)=stifr(n4)+ac
134
135 IF(itherm_fe > 0)THEN
136 ac= fourth*fthe(mc)
137 fthe(n1)=fthe(n1)+ac
138 fthe(n2)=fthe(n2)+ac
139 fthe(n3)=fthe(n3)+ac
140 fthe(n4)=fthe(n4)+ac
141 END IF
142
143 IF(nodadt_therm > 0)THEN
144 ac= fourth*condn(mc)
145 condn(n1)=condn(n1)+ac
146 condn(n2)=condn(n2)+ac
147 condn(n3)=condn(n3)+ac
148 condn(n4)=condn(n4)+ac
149 END IF
150
151 lkinnod(mc)=1
152 stifn(mc)=em20
153 stifr(mc)=em20
154
155 m1=ixc(3,son )
156 IF(lkinnod(m1)==0)THEN
157 lkinnod(m1)=1
158 DO j=1,3
159 a1=half*a(j,m1)
160 a(j,n1)=a(j,n1)+a1
161 a(j,n2)=a(j,n2)+a1
162 END DO
163 a1=half*stifn(m1)
164 stifn(n1)=stifn(n1)+a1
165 stifn(n2)=stifn(n2)+a1
166
167 IF(istatcnd/=0)THEN
168 a1=half*stcont(m1)
169 stcont(n1)=stcont(n1)+a1
170 stcont(n2)=stcont(n2)+a1
171 END IF
172
173 DO j=1,3
174 a1=half*ar(j,m1)
175 ar(j,n1)=ar(j,n1)+a1
176 ar(j,n2)=ar(j,n2)+a1
177 END DO
178 a1=half*stifr(m1)
179 stifr(n1)=stifr(n1)+a1
180 stifr(n2)=stifr(n2)+a1
181
182 IF(itherm_fe > 0)THEN
183 a1= half*fthe(m1)
184 fthe(n1)=fthe(n1)+a1
185 fthe(n2)=fthe(n2)+a1
186 END IF
187
188 IF(nodadt_therm > 0)THEN
189 a1= half*condn(m1)
190 condn(n1)=condn(n1)+a1
191 condn(n2)=condn(n2)+a1
192 END IF
193
194 stifn(m1)=em20
195 stifr(m1)=em20
196 END IF
197
198 m2=ixc(4,son+1)
199 IF(lkinnod(m2)==0)THEN
200 lkinnod(m2)=1
201 DO j=1,3
202 a2=half*a(j,m2)
203 a(j,n2)=a(j,n2)+a2
204 a(j,n3)=a(j,n3)+a2
205 END DO
206 a2=half*stifn(m2)
207 stifn(n2)=stifn(n2)+a2
208 stifn(n3)=stifn(n3)+a2
209
210 IF(istatcnd/=0)THEN
211 a2=half*stcont(m2)
212 stcont(n2)=stcont(n2)+a2
213 stcont(n3)=stcont(n3)+a2
214 END IF
215
216 DO j=1,3
217 a2=half*ar(j,m2)
218 ar(j,n2)=ar(j,n2)+a2
219 ar(j,n3)=ar(j,n3)+a2
220 END DO
221 a2=half*stifr(m2)
222 stifr(n2)=stifr(n2)+a2
223 stifr(n3)=stifr(n3)+a2
224
225 IF(itherm_fe > 0)THEN
226 a2= half*fthe(m2)
227 fthe(n2)=fthe(n2)+a2
228 fthe(n3)=fthe(n3)+a2
229 END IF
230
231 IF(nodadt_therm > 0)THEN
232 a2= half*condn(m2)
233 condn(n2)=condn(n2)+a2
234 condn(n3)=condn(n3)+a2
235 END IF
236
237 stifn(m2)=em20
238 stifr(m2)=em20
239 END IF
240
241 m3=ixc(5,son+2)
242 IF(lkinnod(m3)==0)THEN
243 lkinnod(m3)=1
244 DO j=1,3
245 a3=half*a(j,m3)
246 a(j,n3)=a(j,n3)+a3
247 a(j,n4)=a(j,n4)+a3
248 END DO
249 a3=half*stifn(m3)
250 stifn(n3)=stifn(n3)+a3
251 stifn(n4)=stifn(n4)+a3
252
253 IF(istatcnd/=0)THEN
254 a3=half*stcont(m3)
255 stcont(n3)=stcont(n3)+a3
256 stcont(n4)=stcont(n4)+a3
257 END IF
258
259 DO j=1,3
260 a3=half*ar(j,m3)
261 ar(j,n3)=ar(j,n3)+a3
262 ar(j,n4)=ar(j,n4)+a3
263 END DO
264 a3=half*stifr(m3)
265 stifr(n3)=stifr(n3)+a3
266 stifr(n4)=stifr(n4)+a3
267
268 IF(itherm_fe > 0)THEN
269 a3= half*fthe(m3)
270 fthe(n3)=fthe(n3)+a3
271 fthe(n4)=fthe(n4)+a3
272 END IF
273
274 IF(nodadt_therm > 0)THEN
275 a3= half*condn(m3)
276 condn(n3)=condn(n3)+a3
277 condn(n4)=condn(n4)+a3
278 END IF
279
280 stifn(m3)=em20
281 stifr(m3)=em20
282 END IF
283
284 m4=ixc(2,son+3)
285 IF(lkinnod(m4)==0)THEN
286 lkinnod(m4)=1
287 DO j=1,3
288 a4=half*a(j,m4)
289 a(j,n1)=a(j,n1)+a4
290 a(j,n4)=a(j,n4)+a4
291 END DO
292 a4=half*stifn(m4)
293 stifn(n1)=stifn(n1)+a4
294 stifn(n4)=stifn(n4)+a4
295
296 IF(istatcnd/=0)THEN
297 a4=half*stcont(m4)
298 stcont(n1)=stcont(n1)+a4
299 stcont(n4)=stcont(n4)+a4
300 END IF
301
302 DO j=1,3
303 a4=half*ar(j,m4)
304 ar(j,n1)=ar(j,n1)+a4
305 ar(j,n4)=ar(j,n4)+a4
306 END DO
307 a4=half*stifr(m4)
308 stifr(n1)=stifr(n1)+a4
309 stifr(n4)=stifr(n4)+a4
310
311 IF(itherm_fe > 0)THEN
312 a4= half*fthe(m4)
313 fthe(n1)=fthe(n1)+a4
314 fthe(n4)=fthe(n4)+a4
315 END IF
316
317 IF(nodadt_therm > 0)THEN
318 a4= half*condn(m4)
319 condn(n1)=condn(n1)+a4
320 condn(n4)=condn(n4)+a4
321 END IF
322
323 stifn(m4)=em20
324 stifr(m4)=em20
325 END IF
326
327 END DO
328
329
332
333 son=sh3tree(2,n)
334
335 n1=ixtg(2,n)
336 n2=ixtg(3,n)
337 n3=ixtg(4,n)
338
339 m1=ixtg(4,son+3)
340 IF(lkinnod(m1)==0)THEN
341 lkinnod(m1)=1
342 DO j=1,3
343 a1=half*a(j,m1)
344 a(j,n1)=a(j,n1)+a1
345 a(j,n2)=a(j,n2)+a1
346 END DO
347 a1=half*stifn(m1)
348 stifn(n1)=stifn(n1)+a1
349 stifn(n2)=stifn(n2)+a1
350
351 IF(istatcnd/=0)THEN
352 a1=half*stcont(m1)
353 stcont(n1)=stcont(n1)+a1
354 stcont(n2)=stcont(n2)+a1
355 END IF
356
357 DO j=1,3
358 a1=half*ar(j,m1)
359 ar(j,n1)=ar(j,n1)+a1
360 ar(j,n2)=ar(j,n2)+a1
361 END DO
362 a1=half*stifr(m1)
363 stifr(n1)=stifr(n1)+a1
364 stifr(n2)=stifr(n2)+a1
365
366 IF(itherm_fe > 0)THEN
367 a1= half*fthe(m1)
368 fthe(n1)=fthe(n1)+a1
369 fthe(n2)=fthe(n2)+a1
370 END IF
371
372 IF(nodadt_therm > 0)THEN
373 a1= half*condn(m1)
374 condn(n1)=condn(n1)+a1
375 condn(n2)=condn(n2)+a1
376 END IF
377
378 stifn(m1)=em20
379 stifr(m1)=em20
380 END IF
381
382 m2=ixtg(2,son+3)
383 IF(lkinnod(m2)==0)THEN
384 lkinnod(m2)=1
385 DO j=1,3
386 a2=half*a(j,m2)
387 a(j,n2)=a(j,n2)+a2
388 a(j,n3)=a(j,n3)+a2
389 END DO
390 a2=half*stifn(m2)
391 stifn(n2)=stifn(n2)+a2
392 stifn(n3)=stifn(n3)+a2
393
394 IF(istatcnd/=0)THEN
395 a2=half*stcont(m2)
396 stcont(n2)=stcont(n2)+a2
397 stcont(n3)=stcont(n3)+a2
398 END IF
399
400 DO j=1,3
401 a2=half*ar(j,m2)
402 ar(j,n2)=ar(j,n2)+a2
403 ar(j,n3)=ar(j,n3)+a2
404 END DO
405 a2=half*stifr(m2)
406 stifr(n2)=stifr(n2)+a2
407 stifr(n3)=stifr(n3)+a2
408
409 IF(itherm_fe > 0)THEN
410 a2= half*fthe(m2)
411 fthe(n2)=fthe(n2)+a2
412 fthe(n3)=fthe(n3)+a2
413 END IF
414
415 IF(nodadt_therm > 0)THEN
416 a2= half*condn(m2)
417 condn(n2)=condn(n2)+a2
418 condn(n3)=condn(n3)+a2
419 END IF
420
421 stifn(m2)=em20
422 stifr(m2)=em20
423 END IF
424
425 m3=ixtg(3,son+3)
426 IF(lkinnod(m3)==0)THEN
427 lkinnod(m3)=1
428 DO j=1,3
429 a3=half*a(j,m3)
430 a(j,n3)=a(j,n3)+a3
431 a(j,n1)=a(j,n1)+a3
432 END DO
433 a3=half*stifn(m3)
434 stifn(n3)=stifn(n3)+a3
435 stifn(n1)=stifn(n1)+a3
436
437 IF(istatcnd/=0)THEN
438 a3=half*stcont(m3)
439 stcont(n3)=stcont(n3)+a3
440 stcont(n1)=stcont(n1)+a3
441 END IF
442
443 DO j=1,3
444 a3=half*ar(j,m3)
445 ar(j,n3)=ar(j,n3)+a3
446 ar(j,n1)=ar(j,n1)+a3
447 END DO
448 a3=half*stifr(m3)
449 stifr(n3)=stifr(n3)+a3
450 stifr(n1)=stifr(n1)+a3
451
452 IF(itherm_fe > 0)THEN
453 a3= half*fthe(m3)
454 fthe(n3)=fthe(n3)+a3
455 fthe(n1)=fthe(n1)+a3
456 END IF
457
458 IF(nodadt_therm > 0)THEN
459 a3= half*condn(m3)
460 condn(n3)=condn(n3)+a3
461 condn(n1)=condn(n1)+a3
462 END IF
463
464 stifn(m3)=em20
465 stifr(m3)=em20
466 END IF
467
468 END DO
469
470 END DO
471
472 IF(istatcnd==0) RETURN
473
475
476 acnd(1:3,1:numnod)=a(1:3,1:numnod)
477 arcnd(1:3,1:numnod)=ar(1:3,1:numnod)
478
480 DO nn=1,ll
482
483 n1=ixc(2,n)
484 n2=ixc(3,n)
485 n3=ixc(4,n)
486 n4=ixc(5,n)
487
488
489 rnod(n1)=-one
490 snod(n1)=-one
491 rnod(n2)= one
492 snod(n2)=-one
493 rnod(n3)= one
494 snod(n3)= one
495 rnod(n4)=-one
496 snod(n4)= one
497
498
499 lelt =1
500 nelt(1)=n
501
502 lelt1 =0
503 lelt2 =1
504
505 lev=0
506
507 llnod=0
508 DO WHILE (lev < levelmax)
509 DO le=lelt1+1,lelt2
510
511 ne =nelt(le)
512 IF(sh4tree(3,ne) >= 0) cycle
513
514 m1=ixc(2,ne)
515 m2=ixc(3,ne)
516 m3=ixc(4,ne)
517 m4=ixc(5,ne)
518
519 son=sh4tree(2,ne)
520
521 lelt=lelt+1
522 nelt(lelt)=son
523
524 lelt=lelt+1
525 nelt(lelt)=son+1
526
527 lelt=lelt+1
528 nelt(lelt)=son+2
529
530 lelt=lelt+1
531 nelt(lelt)=son+3
532
533 ni=ixc(3,son)
534 IF(lkinnod(ni)==0.AND.
tagnod(ni)==0)
THEN
535
536
538 llnod=llnod+1
539 lnod(llnod)=ni
540 END IF
541 rnod(ni)=half*(rnod(m1)+rnod(m2))
542 snod(ni)=half*(snod(m1)+snod(m2))
543
544 ni=ixc(4,son+1)
545 IF(lkinnod(ni)==0.AND.
tagnod(ni)==0)
THEN
547 llnod=llnod+1
548 lnod(llnod)=ni
549 END IF
550 rnod(ni)=half*(rnod(m2)+rnod(m3))
551 snod(ni)=half*(snod(m2)+snod(m3))
552
553 ni=ixc(5,son+2)
554 IF(lkinnod(ni)==0.AND.
tagnod(ni)==0)
THEN
556 llnod=llnod+1
557 lnod(llnod)=ni
558 END IF
559 rnod(ni)=half*(rnod(m3)+rnod(m4))
560 snod(ni)=half*(snod(m3)+snod(m4))
561
562 ni=ixc(2,son+3)
563 IF(lkinnod(ni)==0.AND.
tagnod(ni)==0)
THEN
565 llnod=llnod+1
566 lnod(llnod)=ni
567 END IF
568 rnod(ni)=half*(rnod(m4)+rnod(m1))
569 snod(ni)=half*(snod(m4)+snod(m1))
570
571 ni=ixc(4,son)
572 IF(lkinnod(ni)==0)THEN
574 llnod=llnod+1
575 lnod(llnod)=ni
576 END IF
577 rnod(ni)=fourth*(rnod(m1)+rnod(m2)+rnod(m3)+rnod(m4))
578 snod(ni)=fourth*(snod(m1)+snod(m2)+snod(m3)+snod(m4))
579
580 END DO
581
582 lev =lev+1
583 lelt1 =lelt2
584 lelt2 =lelt
585
586 END DO
587
588
589 DO i=1,llnod
590 ni=lnod(i)
591 r =rnod(ni)
592 s =snod(ni)
593 phi =fourth*(one-r)*(one-s)
594 DO j=1,3
595 ac= phi*a(j,ni)
596 a(j,n1)=a(j,n1)+ac
597 END DO
598 stifn(n1)=stifn(n1)+phi*stcont(ni)
599 DO j=1,3
600 ac= phi*ar(j,ni)
601 ar(j,n1)=ar(j,n1)+ac
602 END DO
603 phi=fourth*(one+r)*(one-s)
604 DO j=1,3
605 ac= phi*a(j,ni)
606 a(j,n2)=a(j,n2)+ac
607 END DO
608 stifn(n2)=stifn(n2)+phi*stcont(ni)
609 DO j=1,3
610 ac= phi*ar(j,ni)
611 ar(j,n2)=ar(j,n2)+ac
612 END DO
613 phi=fourth*(one+r)*(one+s)
614 DO j=1,3
615 ac= phi*a(j,ni)
616 a(j,n3)=a(j,n3)+ac
617 END DO
618 stifn(n3)=stifn(n3)+phi*stcont(ni)
619 DO j=1,3
620 ac= phi*ar(j,ni)
621 ar(j,n3)=ar(j,n3)+ac
622 END DO
623 phi=fourth*(one-r)*(one+s)
624 DO j=1,3
625 ac= phi*a(j,ni)
626 a(j,n4)=a(j,n4)+ac
627 END DO
628 stifn(n4)=stifn(n4)+phi*stcont(ni)
629 DO j=1,3
630 ac= phi*ar(j,ni)
631 ar(j,n4)=ar(j,n4)+ac
632 END DO
633 END DO
634
635
636 END DO
637
638
639
641 DO nn=1,ll
643
644 n1=ixtg(2,n)
645 n2=ixtg(3,n)
646 n3=ixtg(4,n)
647
648
649 rnod(n1)= zero
650 snod(n1)= zero
651 rnod(n2)= one
652 snod(n2)= zero
653 rnod(n3)= zero
654 snod(n3)= one
655
656
657 lelt =1
658 nelt(1)=n
659
660 lelt1 =0
661 lelt2 =1
662
663 lev=0
664
665 llnod=0
666 DO WHILE (lev < levelmax)
667 DO le=lelt1+1,lelt2
668
669 ne =nelt(le)
670 IF(sh3tree(3,ne) >= 0) cycle
671
672 m1=ixtg(2,ne)
673 m2=ixtg(3,ne)
674 m3=ixtg(4,ne)
675
676 son=sh3tree(2,ne)
677
678 lelt=lelt+1
679 nelt(lelt)=son
680
681 lelt=lelt+1
682 nelt(lelt)=son+1
683
684 lelt=lelt+1
685 nelt(lelt)=son+2
686
687 lelt=lelt+1
688 nelt(lelt)=son+3
689
690 ni=ixtg(4,son+3)
691 IF(lkinnod(ni)==0.AND.
tagnod(ni)==0)
THEN
693 llnod=llnod+1
694 lnod(llnod)=ni
695 END IF
696 rnod(ni)=half*(rnod(m1)+rnod(m2))
697 snod(ni)=half*(snod(m1)+snod(m2))
698
699 ni=ixtg(2,son+3)
700 IF(lkinnod(ni)==0.AND.
tagnod(ni)==0)
THEN
702 llnod=llnod+1
703 lnod(llnod)=ni
704 END IF
705 rnod(ni)=half*(rnod(m2)+rnod(m3))
706 snod(ni)=half*(snod(m2)+snod(m3))
707
708 ni=ixtg(3,son+3)
709 IF(lkinnod(ni)==0.AND.
tagnod(ni)==0)
THEN
711 llnod=llnod+1
712 lnod(llnod)=ni
713 END IF
714 rnod(ni)=half*(rnod(m3)+rnod(m1))
715 snod(ni)=half*(snod(m3)+snod(m1))
716
717 END DO
718
719 lev =lev+1
720 lelt1 =lelt2
721 lelt2 =lelt
722
723 END DO
724
725
726 DO i=1,llnod
727 ni=lnod(i)
728 r =rnod(ni)
729 s =snod(ni)
730 phi =one-r-s
731 DO j=1,3
732 ac= phi*a(j,ni)
733 a(j,n1)=a(j,n1)+ac
734 END DO
735 stifn(n1)=stifn(n1)+phi*stcont(ni)
736 DO j=1,3
737 ac= phi*ar(j,ni)
738 ar(j,n1)=ar(j,n1)+ac
739 END DO
740 phi=r
741 DO j=1,3
742 ac= phi*a(j,ni)
743 a(j,n2)=a(j,n2)+ac
744 END DO
745 stifn(n2)=stifn(n2)+phi*stcont(ni)
746 DO j=1,3
747 ac= phi*ar(j,ni)
748 ar(j,n2)=ar(j,n2)+ac
749 END DO
750 phi=s
751 DO j=1,3
752 ac= phi*a(j,ni)
753 a(j,n3)=a(j,n3)+ac
754 END DO
755 stifn(n3)=stifn(n3)+phi*stcont(ni)
756 DO j=1,3
757 ac= phi*ar(j,ni)
758 ar(j,n3)=ar(j,n3)+ac
759 END DO
760 END DO
761
762
763 END DO
764
765
766 DEALLOCATE(lnod)
767 DEALLOCATE(nelt)
768 DEALLOCATE(lkinnod)
769 DEALLOCATE(rnod)
770 DEALLOCATE(snod)
771 RETURN
integer, dimension(:), allocatable lsh4upl
integer, dimension(:), allocatable lsh4kin
integer, dimension(:), allocatable lsh3upl
integer, dimension(:), allocatable lsh3kin
integer, dimension(:), allocatable psh4kin
integer, dimension(:), allocatable psh3kin
integer, dimension(:), allocatable psh3upl
integer, dimension(:), allocatable psh4upl
integer, dimension(:), allocatable tagnod