37
38
39
40
41
42 USE elbufdef_mod
43 USE ebcs_mod
45
46
47
48#include "implicit_f.inc"
49
50
51
52#include "param_c.inc"
53#include "com01_c.inc"
54#include "com04_c.inc"
55#include "com08_c.inc"
56#include "scr11_c.inc"
57
58
59
60 INTEGER NSEG,NOD,ISEG(NSEG),LISTE(NOD),IRECT(4,NSEG),
61 . IPARG(NPARG,NGROUP),IELEM(NSEG)
63 . a(3,numnod),x(3,numnod),v(3,numnod),la(3,nod),
64 . p0(nod),vo(nod),po(nod),ms(*),stifn(*),fv(*)
65 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
66 TYPE(t_ebcs_gradp0), INTENT(IN) :: EBCS
67 TYPE(t_segvar) :: SEGVAR
68
69
70
71 INTEGER I,IS,KSEG,ESEG,N1,N2,N3,N4,NG1,NG2,NG3,,,
72 . KLT,KTY,MFT,EAD,,IENER,II(6)
73
75 . orient,rho,c,lcar,roc,alp,fac,
76 . x13,y13,z13,x24,y24,z24,nx,ny,nz,s,
77 . roou,enou,vmx,vmy,vmz,fluxi,fluxo,vn,pn,du,dp,p,
78 . ener,r1,r2,s2,dpdv
79 TYPE(G_BUFEL_) ,POINTER :: GBUF
80
81 irho = ebcs%irho
82 iener = ebcs%iener
83 klt = 0
84 mft = 0
85 IF(irho>0)THEN
86 rho=ebcs%rho*fv(irho)
87 ELSE
88 rho=ebcs%rho
89 ENDIF
90 IF(iener>0)THEN
91 ener=ebcs%ener*fv(iener)
92 ELSE
93 ener=ebcs%ener
94 ENDIF
95 c=ebcs%c
96 lcar=ebcs%lcar
97 r1=ebcs%r1
98 r2=ebcs%r1
99 roc=rho*c
100 alp=zero
101 IF(lcar>0)alp=c*dt1/lcar
102
103
104
105
106 DO i=1,nod
107 la(1,i)=zero
108 la(2,i)=zero
109 la(3,i)=zero
110 ENDDO
111
112 DO is=1,nseg
113 kseg=abs(iseg(is))
114 orient=float(iseg(is)/kseg)
115 eseg=ielem(is)
116 n1=irect(1,is)
117 n2=irect(2,is)
118 n3=irect(3,is)
119 n4=irect(4,is)
120 IF(n4==0 .OR. n4==n3) THEN
121 fac=one_over_6*orient
122 n4=n3
123 ELSE
124 fac=one_over_8*orient
125 ENDIF
126
127 ng1=liste(n1)
128 ng2=liste(n2)
129 ng3=liste(n3)
130 ng4=liste(n4)
131 x13=x(1,ng3)-x(1,ng1)
132 y13=x(2,ng3)-x(2,ng1)
133 z13=x(3,ng3)-x(3,ng1)
134 x24=x(1,ng4)-x(1,ng2)
135 y24=x(2,ng4)-x(2,ng2)
136 z24=x(3,ng4)-x(3,ng2)
137
138 nx=(y13*z24-z13*y24)*fac
139 ny=(z13*x24-x13*z24)*fac
140 nz=(x13*y24-y13*x24)*fac
141
142 la(1,n1)=la(1,n1)+nx
143 la(2,n1)=la(2,n1)+ny
144 la(3,n1)=la(3,n1)+nz
145 la(1,n2)=la(1,n2)+nx
146 la(2,n2)=la(2,n2)+ny
147 la(3,n2)=la(3,n2)+nz
148 la(1,n3)=la(1,n3)+nx
149 la(2,n3)=la(2,n3)+ny
150 la(3,n3)=la(3,n3)+nz
151
152 vmx=v(1,ng1)+v(1,ng2)+v(1,ng3)
153 vmy=v(2,ng1)+v(2,ng2)+v(2,ng3)
154 vmz=v(3,ng1)+v(3,ng2)+v(3,ng3)
155 IF(n4/=n3) THEN
156 la(1,n4)=la(1,n4)+nx
157 la(2,n4)=la(2,n4)+ny
158 la(3,n4)=la(3,n4)+nz
159 vmx=vmx+v(1,ng4)
160 vmy=vmy+v(2,ng4)
161 vmz=vmz+v(3,ng4)
162 ENDIF
163
164
165
166 roou = segvar%RHO(kseg)
167 enou = segvar%EINT(kseg)
168
169 fluxo=(vmx*nx+vmy*ny+vmz*nz)*dt1
170 fluxi=
min(fluxo,zero)
171 fluxo=
max(fluxo,zero)
172
173 dmf=dmf-fluxo*roou-fluxi*rho
174 def=def-fluxo*enou-fluxi*ener
175
176
177
178 segvar%RHO(kseg)=rho
179 segvar%EINT(kseg)=ener
180
181 DO n=1,ngroup
182 kty = iparg(5,n)
183 klt = iparg(2,n)
184 mft = iparg(3,n)
185 IF (kty==1 .AND. eseg<=klt+mft) EXIT
186 ENDDO
187 ead = eseg-mft
188 gbuf => elbuf_tab(n)%GBUF
189 DO i=1,6
190 ii(i) = klt*(i-1)
191 ENDDO
192
193 p =-(gbuf%SIG(ii(1)+ead)+gbuf%SIG(ii(2)+ead)+gbuf%SIG(ii(3)+ead))*third
194
195 p0(n1)=p0(n1)+p*(nx*la(1,n1)+ny*la(2,n1)+nz*la(3,n1))
196 p0(n2)=p0(n2)+p*(nx*la(1,n2)+ny*la(2,n2)+nz*la(3,n2))
197 p0(n3)=p0(n3)+p*(nx*la(1,n3)+ny*la(2,n3)+nz*la(3,n3))
198 IF(n4/=n3) THEN
199 p0(n4)=p0(n4)+p*(nx*la(1,n4)+ny*la(2,n4)+nz*la(3,n4))
200 ENDIF
201
202 ENDDO
203
204
205 DO i=1,nod
206 n=liste(i)
207 s2=la(1,i)**2+la(2,i)**2+la(3,i)**2
208 s=sqrt(s2)
209 vn=(v(1,n)*la(1,i)+v(2,n)*la(2,i)+v(3,n)*la(3,i))/s
210
211 pn=p0(i)/s2+r1*vn+r2*vn*abs(vn)
212 dpdv=roc+r1+two*r2*abs(vn)
213
214 IF(tt>0)THEN
215 du=roc*(vn-vo(i))
216 ELSE
217 du=zero
218 po(i)=pn
219 ENDIF
220 dp=alp*(pn-po(i))
221 p=po(i)+alp*dp+du
222 IF(c==zero)p=pn
223
224
225 a(1,n)=a(1,n)-p*la(1,i)
226 a(2,n)=a(2,n)-p*la(2,i)
227 a(3,n)=a(3,n)-p*la(3,i)
228 stifn(n)=stifn(n)+(two*(s*dpdv)**2)/ms(n)
229
230 def=def-half*(po(i)+p)*dt1*vn*s
231
232 vo(i)=vn
233 po(i)=p
234 ENDDO
235
236 RETURN