34
35
36
38 USE elbufdef_mod
39 USE multi_fvm_mod
40
41
42
43#include "implicit_f.inc"
44
45
46
47#include "com01_c.inc"
48#include "com06_c.inc"
49#include "param_c.inc"
50
51
52
54 INTEGER, INTENT(IN) :: NG
55 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
56 INTEGER, INTENT(IN) :: IPARG(NPARG, *)
57 INTEGER, INTENT(IN) :: ITASK
58 INTEGER, INTENT(IN), TARGET :: IXS(NIXS, *), IXQ(NIXQ, *), IXTG(NIXTG, *)
59 TYPE(MULTI_FVM_STRUCT), INTENT(INOUT) :: MULTI_FVM
60 my_real,
INTENT(IN) :: gravity(4, *)
61 DOUBLE PRECISION,INTENT(INOUT) :: WFEXT
62
63
64
65 TYPE(G_BUFEL_), POINTER :: GBUF
66 INTEGER :: II, I, J, ITY, NEL, NFT, ISOLNOD
67 INTEGER :: IPLA, NODEID, NVERTEX
68 my_real :: rho, etot, vel2, vol, sumflux(5)
69 my_real :: vii(5), gravii(3), pres, vx, vy, vz, wext
70 INTEGER :: NB_FACE, NB_NODE, NODE_LIST(8)
71 INTEGER, DIMENSION(:, :), POINTER :: IX
72 DOUBLE PRECISION :: WFEXTT
73
74
75
76 gbuf =>elbuf_tab(ng)%GBUF
77 nel = iparg(2, ng)
78 nft = iparg(3, ng)
79 ity = iparg(5, ng)
80 isolnod = iparg(28, ng)
81 nb_face = 6
82 wfextt = zero
83 ix =>ixs(1:nixs, 1 + nft : nel + nft)
84 IF (isolnod == 8) THEN
85 nb_node = 8
86 DO j = 1, nb_node
87 node_list = j + 1
88 ENDDO
89 ELSE
90 nb_node = 4
91 node_list(1) = 2
92 node_list(2) = 4
93 node_list(3) = 7
94 node_list(4) = 6
95 ENDIF
96
97 IF (multi_fvm%SYM /= 0) THEN
98 IF (ity == 2) THEN
99
100 nb_face = 4
101 nb_node = 4
102 node_list(1) = 2
103 node_list(2) = 3
104 node_list(3) = 4
105 node_list(4) = 5
106
107 ix => ixq(1:nixq, 1 + nft : nel + nft)
108 ELSEIF (ity == 7) THEN
109
110 nb_face = 3
111 nb_node = 3
112 node_list(1) = 2
113 node_list(2) = 3
114 node_list(3) = 4
115 ix => ixtg(1:nixtg, 1 + nft : nel + nft)
116 ENDIF
117 ENDIF
118
119
120 DO ii = 1, nel
121 i = ii + nft
122
123 vx = gbuf%MOM(ii + 0 * nel)
124 vy = gbuf%MOM(ii + 1 * nel)
125 vz = gbuf%MOM(ii + 2 * nel)
126
127 vel2 = vx * vx + vy * vy + vz * vz
128
129 vii(1) = gbuf%RHO(ii)
130 vii(2) = vii(1) * vx
131 vii(3) = vii(1) * vy
132 vii(4) = vii(1) * vz
133 vii(5) = gbuf%EINT(ii) + half * vii(1) * vel2
134
135
136 vol = gbuf%VOL(ii)
137
138 IF (multi_fvm%SYM == 0 .AND. isolnod /= 4) THEN
139 sumflux(1:5) = multi_fvm%FLUXES(1:5, 1, i) + multi_fvm%FLUXES(1:5, 2, i) +
140 . multi_fvm%FLUXES(1:5, 3, i) + multi_fvm%FLUXES(1:5, 4, i) +
141 . multi_fvm%FLUXES(1:5, 5, i) + multi_fvm%FLUXES(1:5, 6, i)
142 ELSEIF (isolnod == 4) THEN
143 sumflux(1:5) = multi_fvm%FLUXES(1:5, 5, i) + multi_fvm%FLUXES(1:5, 6, i) +
144 . multi_fvm%FLUXES(1:5, 2, i) + multi_fvm%FLUXES(1:5, 4, i)
145 ELSE
146
147 sumflux(1:5) = multi_fvm%FLUXES(1:5, 1, i) + multi_fvm%FLUXES(1:5, 2, i) +
148 . multi_fvm%FLUXES(1:5, 3, i)
149 IF (ity == 2) THEN
150
151 sumflux(1:5) = sumflux(1:5) + multi_fvm%FLUXES(1:5, 4, i)
152 ENDIF
153 ENDIF
154
155 vii(1:5) = vol * vii(1:5) - timestep * sumflux(1:5)
156
157 IF (multi_fvm%SYM == 1) THEN
158 pres = multi_fvm%PRES(i)
159 vii(3) = vii(3) + timestep * gbuf%AREA(ii) * pres
160 ENDIF
161
162 gravii(1:3) = zero
163 nvertex = 0
164
165 DO j = 1, nb_node
166 nodeid = ix(node_list(j), ii)
167 IF(gravity(4, nodeid) == zero) cycle
168 nvertex = nvertex + 1
169 gravii(1) = gravii(1) + gravity(1, nodeid)
170 gravii(2) = gravii(2) + gravity(2, nodeid)
171 gravii(3) = gravii(3) + gravity(3, nodeid)
172 ENDDO
173 IF (nvertex > 0) THEN
174 gravii(1:3) = gravii(1:3) / nvertex
175 ENDIF
176 vii(2:4) = vii(2:4) + timestep * gbuf%RHO(ii) * vol * gravii(1:3)
177 wext = timestep * gbuf%RHO(ii) * vol * (
178 . gravii(1) * multi_fvm%VEL(1, i) +
179 . gravii(2) * multi_fvm%VEL(2, i) +
180 . gravii(3) * multi_fvm%VEL(3, i))
181 vii(5) = vii(5) + wext
182 wfextt = wfextt + wext
183
184 multi_fvm%RHO(i) = vii(1)
185
186 multi_fvm%VEL(1, i) = vii(2)
187 multi_fvm%VEL(2, i) = vii(3)
188 multi_fvm%VEL(3, i) = vii(4)
189
190 multi_fvm%EINT(i) = vii(5)
191 ENDDO
192
193
194
195 wfext=wfext+wfextt
196
197
198
199
200