SPEED
SETUP_DG_ELEM.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
73
74 subroutine setup_dg_elem(nm, sd, tag_mat, cs_nnz_loc, cs_loc, &
75 nn_loc, local_n_num, ne_loc, local_el_num, &
76 xs,ys,zs,&
77 nel_dg_loc, nel_dg_glo, &
78 i4count, mpi_id, mpi_comm, mpi_np,&
79 alfa11, alfa12, alfa13, &
80 alfa21, alfa22, alfa23, &
81 alfa31, alfa32, alfa33, &
82 beta11, beta12, beta13, &
83 beta21, beta22, beta23, &
84 beta31, beta32, beta33, &
85 gamma1, gamma2, gamma3, &
86 delta1, delta2, delta3, &
87 dg_els, scratch_dg_els, &
88 tag_dg_el, tag_dg_yn, tag_dg_link, &
89 tag_dg_frc, val_dg_frc, nload_dg, &
90 con_bc, nface,mpi_file)
91
92
93 use max_var
94 use str_mesh
95 use str_mesh_scratch
96
97
98 implicit none
99
100 include 'SPEED.MPI'
101
102 type(element), dimension(nel_dg_loc), intent(inout) :: dg_els
103 type(scratch_element), dimension(nel_dg_loc), intent(inout) :: scratch_dg_els
104
105 character*70 :: filempi, filename, cmd, mpi_file, filempi_new
106
107 integer*4 :: nm, cs_nnz_loc, nn_loc, ne_loc, nel_dg_loc, nel_dg_glo
108 integer*4 :: mpi_comm, mpi_id, mpi_ierr, mpi_np, nload_dg, tag_ind, nface
109 integer*4 :: im, nn, ie, ned, unitmpi, unitname, nofel, nel_dg_proc
110 integer*4 :: ne1, ne2, ne3, ne4, ic1, ic2, ic3, ic4
111 integer*4 :: ne5, ne6, ne7, ne8, ic5, ic6, ic7, ic8
112 integer*4 :: el_conf, face_conf, face_found, imate, iele, iface
113 integer*4 :: ip, k, j, i, it, ih, ik, tt, indice, node_not_ass
114
115 integer*4, dimension(nm) :: tag_mat, sd
116 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
117 integer*4, dimension(nn_loc) :: local_n_num, i4count
118 integer*4, dimension(ne_loc) :: local_el_num
119 integer*4, dimension(nload_dg) :: tag_dg_el, tag_dg_yn, tag_dg_frc, tag_dg_link
120
121 integer*4, dimension(:), allocatable :: link_glo
122 integer*4, dimension(:,:), allocatable :: ielem_dg
123 integer*4, dimension(:,:), allocatable :: mat_el_face
124 integer*4, dimension(nface,5) :: con_bc
125
126 real*8 :: normal_x, normal_y, normal_z
127 real*8 :: c_alfa11, c_alfa12, c_alfa13, c_alfa21, c_alfa22, c_alfa23, c_alfa31, c_alfa32, c_alfa33
128 real*8 :: c_beta11, c_beta12, c_beta13, c_beta21, c_beta22, c_beta23, c_beta31, c_beta32, c_beta33
129 real*8 :: c_gamma1, c_gamma2, c_gamma3, c_delta1, c_delta2, c_delta3
130 real*8 :: xnod, ynod, znod, csi, eta, zeta, xnod1, ynod1, znod1
131 real*8 :: val1, val2, val3, val4, val5, val6, valmin
132 real*8 :: coef_a, coef_b, coef_t, det_trasf
133
134 real*8, dimension(:), allocatable :: ctgl,wwgl, zn_glo, zt_glo
135 real*8, dimension(:), allocatable :: ct, ww
136 real*8, dimension(nn_loc) :: xs,ys,zs
137 real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13
138 real*8, dimension(ne_loc) :: alfa21,alfa22,alfa23
139 real*8, dimension(ne_loc) :: alfa31,alfa32,alfa33
140 real*8, dimension(ne_loc) :: beta11,beta12,beta13
141 real*8, dimension(ne_loc) :: beta21,beta22,beta23
142 real*8, dimension(ne_loc) :: beta31,beta32,beta33
143 real*8, dimension(ne_loc) :: gamma1,gamma2,gamma3
144 real*8, dimension(ne_loc) :: delta1,delta2,delta3
145 real*8, dimension(nload_dg,2) :: val_dg_frc
146
147 real*8, dimension(:,:), allocatable :: dd
148 real*8, dimension(:,:), allocatable :: normalxyz
149
150 nel_dg_loc = 0
151 ned = cs_loc(0) - 1
152
153
154 do im = 1,nm
155
156 nn = sd(im) +1
157
158 do ie = 1,ned
159 if (cs_loc(cs_loc(ie -1) + 0) .eq. tag_mat(im)) then
160
161 ! face x = - 1
162 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
163 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
164 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
165 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
166
167
168 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
169
170 call get_tag_bc(con_bc, nface, &
171 local_n_num(ne1),local_n_num(ne2), &
172 local_n_num(ne3),local_n_num(ne4), tag_dg_el, nload_dg, tag_ind)
173
174 if (tag_ind == 0) print*, 'ERROR IN GET_TAG_BC'
175
176
177
178 nel_dg_loc = nel_dg_loc + 1
179 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
180 dg_els(nel_dg_loc)%face_el = 1
181 dg_els(nel_dg_loc)%mat = tag_mat(im)
182 dg_els(nel_dg_loc)%spct_deg = nn-1
183 dg_els(nel_dg_loc)%quad_rule = nofqp
184
185 dg_els(nel_dg_loc)%proj_yn = tag_dg_yn(tag_ind)
186 dg_els(nel_dg_loc)%link = tag_dg_link(tag_ind)
187 dg_els(nel_dg_loc)%frac_yn = tag_dg_frc(tag_ind)
188 dg_els(nel_dg_loc)%zt = val_dg_frc(tag_ind,1)
189 dg_els(nel_dg_loc)%zn = val_dg_frc(tag_ind,2)
190
191 call make_normal(1,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
192 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
193
194 dg_els(nel_dg_loc)%nx = normal_x
195 dg_els(nel_dg_loc)%ny = normal_y
196 dg_els(nel_dg_loc)%nz = normal_z
197
198 allocate(ct(dg_els(nel_dg_loc)%quad_rule), ww(dg_els(nel_dg_loc)%quad_rule), &
199 dd(dg_els(nel_dg_loc)%quad_rule, dg_els(nel_dg_loc)%quad_rule))
200
201 call make_lgl_nw(dg_els(nel_dg_loc)%quad_rule, ct, ww, dd)
202
203 allocate(ctgl(dg_els(nel_dg_loc)%quad_rule), wwgl(dg_els(nel_dg_loc)%quad_rule))
204
205 call make_gl_nw(dg_els(nel_dg_loc)%quad_rule, ctgl, wwgl)
206
207 ip = 0
208 do k = 1, dg_els(nel_dg_loc)%quad_rule
209 do j = 1, dg_els(nel_dg_loc)%quad_rule
210 do i = 1,1
211
212 ip = ip + 1
213 scratch_dg_els(nel_dg_loc)%x_nq(ip) = alfa11(ie)*ct(i) + alfa12(ie)*ctgl(j) &
214 + alfa13(ie)*ctgl(k) + beta11(ie)*ctgl(j)*ctgl(k) &
215 + beta12(ie)*ct(i)*ctgl(k) + beta13(ie)*ct(i)*ctgl(j) &
216 + gamma1(ie)*ct(i)*ctgl(j)*ctgl(k) + delta1(ie)
217
218 scratch_dg_els(nel_dg_loc)%y_nq(ip) = alfa21(ie)*ct(i) + alfa22(ie)*ctgl(j) &
219 + alfa23(ie)*ctgl(k) + beta21(ie)*ctgl(j)*ctgl(k) &
220 + beta22(ie)*ct(i)*ctgl(k) + beta23(ie)*ct(i)*ctgl(j) &
221 + gamma2(ie)*ct(i)*ctgl(j)*ctgl(k) + delta2(ie)
222
223 scratch_dg_els(nel_dg_loc)%z_nq(ip) = alfa31(ie)*ct(i) + alfa32(ie)*ctgl(j) &
224 + alfa33(ie)*ctgl(k) + beta31(ie)*ctgl(j)*ctgl(k) &
225 + beta32(ie)*ct(i)*ctgl(k) + beta33(ie)*ct(i)*ctgl(j) &
226 + gamma3(ie)*ct(i)*ctgl(j)*ctgl(k) + delta3(ie)
227
228
229 dg_els(nel_dg_loc)%wx_pl(ip) = 1.d0
230 dg_els(nel_dg_loc)%wy_pl(ip) = wwgl(j)
231 dg_els(nel_dg_loc)%wz_pl(ip) = wwgl(k)
232
233 enddo
234 enddo
235 enddo
236 deallocate(ct,ww,dd)
237 deallocate(ctgl,wwgl)
238
239 endif
240
241 ! face y = - 1
242 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
243 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
244 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
245 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
246
247
248 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
249
250
251 call get_tag_bc(con_bc, nface, &
252 local_n_num(ne1),local_n_num(ne2), &
253 local_n_num(ne3),local_n_num(ne4), tag_dg_el, nload_dg, tag_ind)
254
255 if (tag_ind == 0) print*, 'ERROR IN GET_TAG_BC'
256
257
258 nel_dg_loc = nel_dg_loc +1
259 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
260 dg_els(nel_dg_loc)%face_el = 2
261 dg_els(nel_dg_loc)%mat = tag_mat(im)
262 dg_els(nel_dg_loc)%spct_deg = nn-1
263 dg_els(nel_dg_loc)%quad_rule = nofqp
264
265 dg_els(nel_dg_loc)%proj_yn = tag_dg_yn(tag_ind)
266 dg_els(nel_dg_loc)%link = tag_dg_link(tag_ind)
267 dg_els(nel_dg_loc)%frac_yn = tag_dg_frc(tag_ind)
268 dg_els(nel_dg_loc)%zt = val_dg_frc(tag_ind,1)
269 dg_els(nel_dg_loc)%zn = val_dg_frc(tag_ind,2)
270
271
272 call make_normal(2,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
273 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
274
275 dg_els(nel_dg_loc)%nx = normal_x
276 dg_els(nel_dg_loc)%ny = normal_y
277 dg_els(nel_dg_loc)%nz = normal_z
278
279 allocate(ct(dg_els(nel_dg_loc)%quad_rule), ww(dg_els(nel_dg_loc)%quad_rule), &
280 dd(dg_els(nel_dg_loc)%quad_rule,dg_els(nel_dg_loc)%quad_rule))
281
282 call make_lgl_nw(dg_els(nel_dg_loc)%quad_rule,ct,ww,dd)
283
284 allocate(ctgl(dg_els(nel_dg_loc)%quad_rule),wwgl(dg_els(nel_dg_loc)%quad_rule))
285
286 call make_gl_nw(dg_els(nel_dg_loc)%quad_rule,ctgl,wwgl)
287
288 ip = 0
289 do k = 1, dg_els(nel_dg_loc)%quad_rule
290 do j = 1,1
291 do i = 1, dg_els(nel_dg_loc)%quad_rule
292 ip = ip + 1
293 scratch_dg_els(nel_dg_loc)%x_nq(ip) = alfa11(ie)*ctgl(i) + alfa12(ie)*ct(j) &
294 + alfa13(ie)*ctgl(k) + beta11(ie)*ct(j)*ctgl(k) &
295 + beta12(ie)*ctgl(i)*ctgl(k) + beta13(ie)*ctgl(i)*ct(j) &
296 + gamma1(ie)*ctgl(i)*ct(j)*ctgl(k) + delta1(ie)
297
298 scratch_dg_els(nel_dg_loc)%y_nq(ip) = alfa21(ie)*ctgl(i) + alfa22(ie)*ct(j) &
299 + alfa23(ie)*ctgl(k) + beta21(ie)*ct(j)*ctgl(k) &
300 + beta22(ie)*ctgl(i)*ctgl(k) + beta23(ie)*ctgl(i)*ct(j) &
301 + gamma2(ie)*ctgl(i)*ct(j)*ctgl(k) + delta2(ie)
302
303 scratch_dg_els(nel_dg_loc)%z_nq(ip) = alfa31(ie)*ctgl(i) + alfa32(ie)*ct(j) &
304 + alfa33(ie)*ctgl(k) + beta31(ie)*ct(j)*ctgl(k) &
305 + beta32(ie)*ctgl(i)*ctgl(k) + beta33(ie)*ctgl(i)*ct(j) &
306 + gamma3(ie)*ctgl(i)*ct(j)*ctgl(k) + delta3(ie)
307
308 dg_els(nel_dg_loc)%wx_pl(ip) = wwgl(i)
309 dg_els(nel_dg_loc)%wy_pl(ip) = 1.d0
310 dg_els(nel_dg_loc)%wz_pl(ip) = wwgl(k)
311
312 enddo
313 enddo
314 enddo
315 deallocate(ct,ww,dd)
316 deallocate(ctgl,wwgl)
317 endif
318
319 ! face z = - 1
320 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
321 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
322 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
323 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
324
325
326 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
327
328 call get_tag_bc(con_bc, nface, &
329 local_n_num(ne1),local_n_num(ne2), &
330 local_n_num(ne3),local_n_num(ne4), tag_dg_el, nload_dg, tag_ind)
331
332 if (tag_ind == 0) print*, 'ERROR IN GET_TAG_BC'
333
334
335 nel_dg_loc = nel_dg_loc +1
336 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
337 dg_els(nel_dg_loc)%face_el = 3
338 dg_els(nel_dg_loc)%mat = tag_mat(im)
339 dg_els(nel_dg_loc)%spct_deg = nn-1
340 dg_els(nel_dg_loc)%quad_rule = nofqp
341
342 dg_els(nel_dg_loc)%proj_yn = tag_dg_yn(tag_ind)
343 dg_els(nel_dg_loc)%link = tag_dg_link(tag_ind)
344 dg_els(nel_dg_loc)%frac_yn = tag_dg_frc(tag_ind)
345 dg_els(nel_dg_loc)%zt = val_dg_frc(tag_ind,1)
346 dg_els(nel_dg_loc)%zn = val_dg_frc(tag_ind,2)
347
348 call make_normal(3,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
349 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
350
351 dg_els(nel_dg_loc)%nx = normal_x
352 dg_els(nel_dg_loc)%ny = normal_y
353 dg_els(nel_dg_loc)%nz = normal_z
354
355 allocate(ct(dg_els(nel_dg_loc)%quad_rule), ww(dg_els(nel_dg_loc)%quad_rule), &
356 dd(dg_els(nel_dg_loc)%quad_rule,dg_els(nel_dg_loc)%quad_rule))
357
358 call make_lgl_nw(dg_els(nel_dg_loc)%quad_rule, ct, ww, dd)
359
360 allocate(ctgl(dg_els(nel_dg_loc)%quad_rule), wwgl(dg_els(nel_dg_loc)%quad_rule))
361
362 call make_gl_nw(dg_els(nel_dg_loc)%quad_rule, ctgl, wwgl)
363
364 ip = 0
365 do k = 1,1
366 do j = 1, dg_els(nel_dg_loc)%quad_rule
367 do i = 1, dg_els(nel_dg_loc)%quad_rule
368 ip = ip + 1
369 scratch_dg_els(nel_dg_loc)%x_nq(ip) = alfa11(ie)*ctgl(i) + alfa12(ie)*ctgl(j) &
370 + alfa13(ie)*ct(k) + beta11(ie)*ctgl(j)*ct(k) &
371 + beta12(ie)*ctgl(i)*ct(k) + beta13(ie)*ctgl(i)*ctgl(j) &
372 + gamma1(ie)*ctgl(i)*ctgl(j)*ct(k) + delta1(ie)
373
374 scratch_dg_els(nel_dg_loc)%y_nq(ip) = alfa21(ie)*ctgl(i) + alfa22(ie)*ctgl(j) &
375 + alfa23(ie)*ct(k) + beta21(ie)*ctgl(j)*ct(k) &
376 + beta22(ie)*ctgl(i)*ct(k) + beta23(ie)*ctgl(i)*ctgl(j) &
377 + gamma2(ie)*ctgl(i)*ctgl(j)*ct(k) + delta2(ie)
378
379 scratch_dg_els(nel_dg_loc)%z_nq(ip) = alfa31(ie)*ctgl(i) + alfa32(ie)*ctgl(j) &
380 + alfa33(ie)*ct(k) + beta31(ie)*ctgl(j)*ct(k) &
381 + beta32(ie)*ctgl(i)*ct(k) + beta33(ie)*ctgl(i)*ctgl(j) &
382 + gamma3(ie)*ctgl(i)*ctgl(j)*ct(k) + delta3(ie)
383
384
385 dg_els(nel_dg_loc)%wx_pl(ip) = wwgl(i)
386 dg_els(nel_dg_loc)%wy_pl(ip) = wwgl(j)
387 dg_els(nel_dg_loc)%wz_pl(ip) = 1.d0
388
389 enddo
390 enddo
391 enddo
392 deallocate(ct,ww,dd)
393 deallocate(ctgl,wwgl)
394 endif
395
396
397 ! face x = 1
398 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
399 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
400 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
401 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
402
403
404 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
405
406 call get_tag_bc(con_bc, nface, &
407 local_n_num(ne1),local_n_num(ne2), &
408 local_n_num(ne3),local_n_num(ne4), tag_dg_el, nload_dg, tag_ind)
409
410 if (tag_ind == 0) print*, 'ERROR IN GET_TAG_BC'
411
412 nel_dg_loc = nel_dg_loc +1
413 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
414 dg_els(nel_dg_loc)%face_el = 4
415 dg_els(nel_dg_loc)%mat = tag_mat(im)
416 dg_els(nel_dg_loc)%spct_deg = nn-1
417 dg_els(nel_dg_loc)%quad_rule = nofqp
418
419 dg_els(nel_dg_loc)%proj_yn = tag_dg_yn(tag_ind)
420 dg_els(nel_dg_loc)%link = tag_dg_link(tag_ind)
421 dg_els(nel_dg_loc)%frac_yn = tag_dg_frc(tag_ind)
422 dg_els(nel_dg_loc)%zt = val_dg_frc(tag_ind,1)
423 dg_els(nel_dg_loc)%zn = val_dg_frc(tag_ind,2)
424
425 call make_normal(4,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
426 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
427
428 dg_els(nel_dg_loc)%nx = normal_x
429 dg_els(nel_dg_loc)%ny = normal_y
430 dg_els(nel_dg_loc)%nz = normal_z
431
432 allocate(ct(dg_els(nel_dg_loc)%quad_rule),ww(dg_els(nel_dg_loc)%quad_rule), &
433 dd(dg_els(nel_dg_loc)%quad_rule,dg_els(nel_dg_loc)%quad_rule))
434
435 call make_lgl_nw(dg_els(nel_dg_loc)%quad_rule, ct, ww, dd)
436
437 allocate(ctgl(dg_els(nel_dg_loc)%quad_rule),wwgl(dg_els(nel_dg_loc)%quad_rule))
438
439 call make_gl_nw(dg_els(nel_dg_loc)%quad_rule, ctgl, wwgl)
440
441 ip = 0
442 do k = 1, dg_els(nel_dg_loc)%quad_rule
443 do j = 1, dg_els(nel_dg_loc)%quad_rule
444 do i = dg_els(nel_dg_loc)%quad_rule, dg_els(nel_dg_loc)%quad_rule
445 ip = ip + 1
446 scratch_dg_els(nel_dg_loc)%x_nq(ip) = alfa11(ie)*ct(i) + alfa12(ie)*ctgl(j) &
447 + alfa13(ie)*ctgl(k) + beta11(ie)*ctgl(j)*ctgl(k) &
448 + beta12(ie)*ct(i)*ctgl(k) + beta13(ie)*ct(i)*ctgl(j) &
449 + gamma1(ie)*ct(i)*ctgl(j)*ctgl(k) + delta1(ie)
450
451 scratch_dg_els(nel_dg_loc)%y_nq(ip) = alfa21(ie)*ct(i) + alfa22(ie)*ctgl(j) &
452 + alfa23(ie)*ctgl(k) + beta21(ie)*ctgl(j)*ctgl(k) &
453 + beta22(ie)*ct(i)*ctgl(k) + beta23(ie)*ct(i)*ctgl(j) &
454 + gamma2(ie)*ct(i)*ctgl(j)*ctgl(k) + delta2(ie)
455
456 scratch_dg_els(nel_dg_loc)%z_nq(ip) = alfa31(ie)*ct(i) + alfa32(ie)*ctgl(j) &
457 + alfa33(ie)*ctgl(k) + beta31(ie)*ctgl(j)*ctgl(k) &
458 + beta32(ie)*ct(i)*ctgl(k) + beta33(ie)*ct(i)*ctgl(j) &
459 + gamma3(ie)*ct(i)*ctgl(j)*ctgl(k) + delta3(ie)
460
461 dg_els(nel_dg_loc)%wx_pl(ip) = 1.d0
462 dg_els(nel_dg_loc)%wy_pl(ip) = wwgl(j)
463 dg_els(nel_dg_loc)%wz_pl(ip) = wwgl(k)
464
465 enddo
466 enddo
467 enddo
468 deallocate(ct,ww,dd)
469 deallocate(ctgl,wwgl)
470 endif
471
472
473 ! face y = 1
474 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
475 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
476 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
477 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
478
479
480 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
481
482 call get_tag_bc(con_bc, nface, &
483 local_n_num(ne1),local_n_num(ne2), &
484 local_n_num(ne3),local_n_num(ne4), tag_dg_el, nload_dg, tag_ind)
485
486 if (tag_ind == 0) print*, 'ERROR IN GET_TAG_BC'
487
488 nel_dg_loc = nel_dg_loc +1
489 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
490 dg_els(nel_dg_loc)%face_el = 5
491 dg_els(nel_dg_loc)%mat = tag_mat(im)
492 dg_els(nel_dg_loc)%spct_deg = nn-1
493 dg_els(nel_dg_loc)%quad_rule = nofqp
494
495 dg_els(nel_dg_loc)%proj_yn = tag_dg_yn(tag_ind)
496 dg_els(nel_dg_loc)%link = tag_dg_link(tag_ind)
497 dg_els(nel_dg_loc)%frac_yn = tag_dg_frc(tag_ind)
498 dg_els(nel_dg_loc)%zt = val_dg_frc(tag_ind,1)
499 dg_els(nel_dg_loc)%zn = val_dg_frc(tag_ind,2)
500
501 call make_normal(5,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
502 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
503
504 dg_els(nel_dg_loc)%nx = normal_x
505 dg_els(nel_dg_loc)%ny = normal_y
506 dg_els(nel_dg_loc)%nz = normal_z
507
508
509 allocate(ct(dg_els(nel_dg_loc)%quad_rule), ww(dg_els(nel_dg_loc)%quad_rule), &
510 dd(dg_els(nel_dg_loc)%quad_rule,dg_els(nel_dg_loc)%quad_rule))
511
512 call make_lgl_nw(dg_els(nel_dg_loc)%quad_rule, ct, ww, dd)
513
514 allocate(ctgl(dg_els(nel_dg_loc)%quad_rule),wwgl(dg_els(nel_dg_loc)%quad_rule))
515
516 call make_gl_nw(dg_els(nel_dg_loc)%quad_rule, ctgl, wwgl)
517 ip = 0
518 do k = 1, dg_els(nel_dg_loc)%quad_rule
519 do j = dg_els(nel_dg_loc)%quad_rule, dg_els(nel_dg_loc)%quad_rule
520 do i = 1, dg_els(nel_dg_loc)%quad_rule
521 ip = ip + 1
522 scratch_dg_els(nel_dg_loc)%x_nq(ip) = alfa11(ie)*ctgl(i) + alfa12(ie)*ct(j) &
523 + alfa13(ie)*ctgl(k) + beta11(ie)*ct(j)*ctgl(k) &
524 + beta12(ie)*ctgl(i)*ctgl(k) + beta13(ie)*ctgl(i)*ct(j) &
525 + gamma1(ie)*ctgl(i)*ct(j)*ctgl(k) + delta1(ie)
526
527 scratch_dg_els(nel_dg_loc)%y_nq(ip) = alfa21(ie)*ctgl(i) + alfa22(ie)*ct(j) &
528 + alfa23(ie)*ctgl(k) + beta21(ie)*ct(j)*ctgl(k) &
529 + beta22(ie)*ctgl(i)*ctgl(k) + beta23(ie)*ctgl(i)*ct(j) &
530 + gamma2(ie)*ctgl(i)*ct(j)*ctgl(k) + delta2(ie)
531
532 scratch_dg_els(nel_dg_loc)%z_nq(ip) = alfa31(ie)*ctgl(i) + alfa32(ie)*ct(j) &
533 + alfa33(ie)*ctgl(k) + beta31(ie)*ct(j)*ctgl(k) &
534 + beta32(ie)*ctgl(i)*ctgl(k) + beta33(ie)*ctgl(i)*ct(j) &
535 + gamma3(ie)*ctgl(i)*ct(j)*ctgl(k) + delta3(ie)
536
537 dg_els(nel_dg_loc)%wx_pl(ip) = wwgl(i)
538 dg_els(nel_dg_loc)%wy_pl(ip) = 1.d0
539 dg_els(nel_dg_loc)%wz_pl(ip) = wwgl(k)
540
541 enddo
542 enddo
543 enddo
544 deallocate(ct,ww,dd)
545 deallocate(ctgl,wwgl)
546 endif
547
548
549 ! face z = 1
550 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
551 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
552 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
553 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
554
555
556 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
557
558 call get_tag_bc(con_bc, nface, &
559 local_n_num(ne1),local_n_num(ne2), &
560 local_n_num(ne3),local_n_num(ne4), tag_dg_el, nload_dg, tag_ind)
561
562 if (tag_ind == 0) print*, 'ERROR IN GET_TAG_BC'
563
564 nel_dg_loc = nel_dg_loc + 1
565 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
566 dg_els(nel_dg_loc)%face_el = 6
567 dg_els(nel_dg_loc)%mat = tag_mat(im)
568 dg_els(nel_dg_loc)%spct_deg = nn-1
569 dg_els(nel_dg_loc)%quad_rule = nofqp
570
571 dg_els(nel_dg_loc)%proj_yn = tag_dg_yn(tag_ind)
572 dg_els(nel_dg_loc)%link = tag_dg_link(tag_ind)
573 dg_els(nel_dg_loc)%frac_yn = tag_dg_frc(tag_ind)
574 dg_els(nel_dg_loc)%zt = val_dg_frc(tag_ind,1)
575 dg_els(nel_dg_loc)%zn = val_dg_frc(tag_ind,2)
576
577 call make_normal(6,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
578 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
579
580 dg_els(nel_dg_loc)%nx = normal_x
581 dg_els(nel_dg_loc)%ny = normal_y
582 dg_els(nel_dg_loc)%nz = normal_z
583
584 allocate(ct(dg_els(nel_dg_loc)%quad_rule),ww(dg_els(nel_dg_loc)%quad_rule), &
585 dd(dg_els(nel_dg_loc)%quad_rule,dg_els(nel_dg_loc)%quad_rule))
586
587 call make_lgl_nw(dg_els(nel_dg_loc)%quad_rule, ct, ww, dd)
588
589 allocate(ctgl(dg_els(nel_dg_loc)%quad_rule), wwgl(dg_els(nel_dg_loc)%quad_rule))
590
591 call make_gl_nw(dg_els(nel_dg_loc)%quad_rule, ctgl, wwgl)
592
593 ip = 0
594 do k = dg_els(nel_dg_loc)%quad_rule,dg_els(nel_dg_loc)%quad_rule
595 do j = 1,dg_els(nel_dg_loc)%quad_rule
596 do i = 1,dg_els(nel_dg_loc)%quad_rule
597 ip = ip + 1
598 scratch_dg_els(nel_dg_loc)%x_nq(ip) = alfa11(ie)*ctgl(i) + alfa12(ie)*ctgl(j) &
599 + alfa13(ie)*ct(k) + beta11(ie)*ctgl(j)*ct(k) &
600 + beta12(ie)*ctgl(i)*ct(k) + beta13(ie)*ctgl(i)*ctgl(j) &
601 + gamma1(ie)*ctgl(i)*ctgl(j)*ct(k) + delta1(ie)
602
603 scratch_dg_els(nel_dg_loc)%y_nq(ip) = alfa21(ie)*ctgl(i) + alfa22(ie)*ctgl(j) &
604 + alfa23(ie)*ct(k) + beta21(ie)*ctgl(j)*ct(k) &
605 + beta22(ie)*ctgl(i)*ct(k) + beta23(ie)*ctgl(i)*ctgl(j) &
606 + gamma2(ie)*ctgl(i)*ctgl(j)*ct(k) + delta2(ie)
607
608 scratch_dg_els(nel_dg_loc)%z_nq(ip) = alfa31(ie)*ctgl(i) + alfa32(ie)*ctgl(j) &
609 + alfa33(ie)*ct(k) + beta31(ie)*ctgl(j)*ct(k) &
610 + beta32(ie)*ctgl(i)*ct(k) + beta33(ie)*ctgl(i)*ctgl(j) &
611 + gamma3(ie)*ctgl(i)*ctgl(j)*ct(k) + delta3(ie)
612
613 dg_els(nel_dg_loc)%wx_pl(ip) = wwgl(i)
614 dg_els(nel_dg_loc)%wy_pl(ip) = wwgl(j)
615 dg_els(nel_dg_loc)%wz_pl(ip) = 1.d0
616
617 enddo
618 enddo
619 enddo
620 deallocate(ct,ww,dd)
621 deallocate(ctgl,wwgl)
622 endif
623 endif
624 enddo
625 enddo
626
627
628 !write in order
629 !- material - element - face - normalx, normaly, normalz
630
631 filempi = 'NORM000000.mpi'
632 unitmpi = 40
633
634 if (mpi_id .lt. 10) then
635 write(filempi(10:10),'(i1)') mpi_id
636 else if (mpi_id .lt. 100) then
637 write(filempi(9:10),'(i2)') mpi_id
638 else if (mpi_id .lt. 1000) then
639 write(filempi(8:10),'(i3)') mpi_id
640 else if (mpi_id .lt. 10000) then
641 write(filempi(7:10),'(i4)') mpi_id
642 else if (mpi_id .lt. 100000) then
643 write(filempi(6:10),'(i5)') mpi_id
644 else if (mpi_id .lt. 1000000) then
645 write(filempi(5:10),'(i6)') mpi_id
646 endif
647
648 if(len_trim(mpi_file) .ne. 70) then
649 filempi_new = mpi_file(1:len_trim(mpi_file)) // '/' // filempi
650 else
651 filempi_new = filempi
652 endif
653
654
655 open(unitmpi,file=filempi_new)
656 write(unitmpi,*) nel_dg_loc
657 do i = 1, nel_dg_loc
658 write(unitmpi,"(1I2,1X,1I12,1X,1I2,1X,3(1X,ES12.4),1X,1I2,2(1X,ES12.4),1X,1I3)") &
659 dg_els(i)%mat, dg_els(i)%ind_el, dg_els(i)%face_el, &
660 dg_els(i)%nx, dg_els(i)%ny, dg_els(i)%nz, dg_els(i)%frac_yn, dg_els(i)%zt, dg_els(i)%zn, dg_els(i)%link
661 enddo
662 close(unitmpi)
663
664 call mpi_barrier(mpi_comm, mpi_ierr)
665
666
667
668 if(mpi_id .eq. 0) then
669
670 filename = 'NORMALL.input'
671 unitname = 500
672 open(unitname,file=filename)
673 write(unitname,*) nel_dg_glo
674 !close(unitname)
675
676 allocate(mat_el_face(nel_dg_glo,4), normalxyz(nel_dg_glo,3), zt_glo(nel_dg_glo), zn_glo(nel_dg_glo), &
677 link_glo(nel_dg_glo))
678 k = 1
679
680 do i = 1, mpi_np
681
682 filempi = 'NORM000000.mpi'
683
684 if (i-1 .lt. 10) then
685 write(filempi(10:10),'(i1)') i-1
686 else if (i-1 .lt. 100) then
687 write(filempi(9:10),'(i2)') i-1
688 else if (i-1 .lt. 1000) then
689 write(filempi(8:10),'(i3)') i-1
690 else if (i-1 .lt. 10000) then
691 write(filempi(7:10),'(i4)') i-1
692 else if (i-1 .lt. 100000) then
693 write(filempi(6:10),'(i5)') i-1
694 else if (i-1 .lt. 1000000) then
695 write(filempi(5:10),'(i6)') i-1
696 endif
697
698 if(len_trim(mpi_file) .ne. 70) then
699 filempi_new = mpi_file(1:len_trim(mpi_file)) // '/' // filempi
700 else
701 filempi_new = filempi
702 endif
703
704 open(unitmpi,file=filempi_new)
705 read(unitmpi,*) nel_dg_proc
706
707 do j = 1, nel_dg_proc
708 read(unitmpi,*) mat_el_face(k,1),mat_el_face(k,2), mat_el_face(k,3), &
709 normalxyz(k,1), normalxyz(k,2), normalxyz(k,3), mat_el_face(k,4), zt_glo(k), zn_glo(k), link_glo(k)
710 k=k+1
711 enddo
712
713 close(unitmpi)
714
715 enddo
716
717 do j = 1, nel_dg_glo
718 write(unitname,*) mat_el_face(j,1),mat_el_face(j,2), mat_el_face(j,3), &
719 normalxyz(j,1), normalxyz(j,2), normalxyz(j,3), mat_el_face(j,4), zt_glo(j), zn_glo(j), link_glo(j)
720 enddo
721
722 deallocate(mat_el_face, normalxyz, zn_glo, zt_glo)
723 close(unitname)
724
725
726 endif
727
728
729 call mpi_barrier(mpi_comm, mpi_ierr)
730
731
732
733
734
735 return
736
737 end subroutine setup_dg_elem
get_tag_bc
subroutine get_tag_bc(cn_bc, nfac, v1, v2, v3, v4, tag_dg, nl_dg, ind)
Findes index of a boundary element.
Definition GET_TAG_BC.f90:34
make_gl_nw
subroutine make_gl_nw(ngp, xabsc, weig)
Makes Gauss Legendre nodes and weights.
Definition MAKE_GL_NW.f90:29
make_lgl_nw
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
Definition MAKE_LGL_NW.f90:29
make_normal
subroutine make_normal(ind, xs1, xs2, xs3, xs4, ys1, ys2, ys3, ys4, zs1, zs2, zs3, zs4, nx, ny, nz, par, arr)
Makes normal vector of a given surface.
Definition MAKE_NORMAL.f90:48
setup_dg_elem
subroutine setup_dg_elem(nm, sd, tag_mat, cs_nnz_loc, cs_loc, nn_loc, local_n_num, ne_loc, local_el_num, xs, ys, zs, nel_dg_loc, nel_dg_glo, i4count, mpi_id, mpi_comm, mpi_np, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, dg_els, scratch_dg_els, tag_dg_el, tag_dg_yn, tag_dg_link, tag_dg_frc, val_dg_frc, nload_dg, con_bc, nface, mpi_file)
Setup for DG data structure.
Definition SETUP_DG_ELEM.f90:91
max_var
Set maximal bounds.
Definition MODULES.f90:54
max_var::nofqp
integer, parameter nofqp
max number of 1-D quadrature point per element
Definition MODULES.f90:57
str_mesh_scratch
Contains mesh structure (scratch)
Definition MODULES.f90:67
str_mesh
Contains mesh structure for DG interface elements.
Definition MODULES.f90:84
str_mesh::element
Definition MODULES.f90:88
str_mesh_scratch::scratch_element
Definition MODULES.f90:71