SPEED
SETUP_DG.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
68
69 subroutine setup_dg(nm, sd, tag_mat, cs_nnz_loc, cs_loc, &
70 nn_loc, local_n_num, ne_loc, local_el_num, &
71 xs,ys,zs,&
72 nel_dg_loc, nel_dg_glo, &
73 i4count, mpi_id, mpi_comm, &
74 alfa11, alfa12, alfa13, &
75 alfa21, alfa22, alfa23, &
76 alfa31, alfa32, alfa33, &
77 beta11, beta12, beta13, &
78 beta21, beta22, beta23, &
79 beta31, beta32, beta33, &
80 gamma1, gamma2, gamma3, &
81 delta1, delta2, delta3, &
82 faces, area_nodes)
83
84
85 implicit none
86
87 include 'SPEED.MPI'
88
89 integer*4 :: nm, cs_nnz_loc, nn_loc, ne_loc, nel_dg_loc, nel_dg_glo
90 integer*4 :: im, nn, ie, ned, err_out
91 integer*4 :: ne1, ne2, ne3, ne4, ic1, ic2, ic3, ic4
92 integer*4 :: ne5, ne6, ne7, ne8, ic5, ic6, ic7, ic8
93 integer*4 :: mpi_comm, mpi_id, mpi_ierr
94
95 integer*4, dimension(nm) :: tag_mat, sd
96 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
97 integer*4, dimension(nn_loc) :: local_n_num, i4count
98 integer*4, dimension(ne_loc) :: local_el_num
99
100 integer*4, dimension(3,nel_dg_loc), intent(inout) :: faces
101
102 real*8 :: surf
103
104 real*8, dimension(nn_loc) :: xs,ys,zs
105 real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13
106 real*8, dimension(ne_loc) :: alfa21,alfa22,alfa23
107 real*8, dimension(ne_loc) :: alfa31,alfa32,alfa33
108 real*8, dimension(ne_loc) :: beta11,beta12,beta13
109 real*8, dimension(ne_loc) :: beta21,beta22,beta23
110 real*8, dimension(ne_loc) :: beta31,beta32,beta33
111 real*8, dimension(ne_loc) :: gamma1,gamma2,gamma3
112 real*8, dimension(ne_loc) :: delta1,delta2,delta3
113
114 real*8, dimension(25,nel_dg_loc), intent(inout) :: area_nodes
115
116
117!*****************************************************************************************
118!loading of data structure faces containing the description for DG faces
119!*****************************************************************************************
120
121 nel_dg_loc = 0
122 ned = cs_loc(0) - 1
123
124 do im = 1,nm
125
126 nn = sd(im) +1
127
128 do ie = 1,ned
129 if (cs_loc(cs_loc(ie -1) + 0) .eq. tag_mat(im)) then
130
131 ! face x = - 1
132 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
133 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
134 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
135 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
136
137
138 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
139 nel_dg_loc = nel_dg_loc +1
140
141 call get_area_face(xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
142 zs(ne1), zs(ne2), zs(ne3), zs(ne4), surf, ie)
143
144 ! face x = 1
145
146 faces(1,nel_dg_loc) = tag_mat(im)
147 faces(2,nel_dg_loc) = local_el_num(ie)
148 faces(3,nel_dg_loc) = 1
149
150 area_nodes(1,nel_dg_loc) = surf
151 area_nodes(2,nel_dg_loc) = alfa11(ie)
152 area_nodes(3,nel_dg_loc) = alfa12(ie)
153 area_nodes(4,nel_dg_loc) = alfa13(ie)
154
155 area_nodes(5,nel_dg_loc) = alfa21(ie)
156 area_nodes(6,nel_dg_loc) = alfa22(ie)
157 area_nodes(7,nel_dg_loc) = alfa23(ie)
158
159 area_nodes(8,nel_dg_loc) = alfa31(ie)
160 area_nodes(9,nel_dg_loc) = alfa32(ie)
161 area_nodes(10,nel_dg_loc) = alfa33(ie)
162
163 area_nodes(11,nel_dg_loc) = beta11(ie)
164 area_nodes(12,nel_dg_loc) = beta12(ie)
165 area_nodes(13,nel_dg_loc) = beta13(ie)
166
167 area_nodes(14,nel_dg_loc) = beta21(ie)
168 area_nodes(15,nel_dg_loc) = beta22(ie)
169 area_nodes(16,nel_dg_loc) = beta23(ie)
170
171 area_nodes(17,nel_dg_loc) = beta31(ie)
172 area_nodes(18,nel_dg_loc) = beta32(ie)
173 area_nodes(19,nel_dg_loc) = beta33(ie)
174
175 area_nodes(20,nel_dg_loc) = gamma1(ie)
176 area_nodes(21,nel_dg_loc) = gamma2(ie)
177 area_nodes(22,nel_dg_loc) = gamma3(ie)
178
179 area_nodes(23,nel_dg_loc) = delta1(ie)
180 area_nodes(24,nel_dg_loc) = delta2(ie)
181 area_nodes(25,nel_dg_loc) = delta3(ie)
182
183 endif
184
185 ! face y = - 1
186 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
187 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
188 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
189 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
190
191
192 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
193 nel_dg_loc = nel_dg_loc +1
194
195 call get_area_face(xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
196 zs(ne1), zs(ne2), zs(ne3), zs(ne4), surf, ie)
197
198 ! face y = 1
199
200 faces(1,nel_dg_loc) = tag_mat(im)
201 faces(2,nel_dg_loc) = local_el_num(ie)
202 faces(3,nel_dg_loc) = 2
203
204 area_nodes(1,nel_dg_loc) = surf
205 area_nodes(2,nel_dg_loc) = alfa11(ie)
206 area_nodes(3,nel_dg_loc) = alfa12(ie)
207 area_nodes(4,nel_dg_loc) = alfa13(ie)
208
209 area_nodes(5,nel_dg_loc) = alfa21(ie)
210 area_nodes(6,nel_dg_loc) = alfa22(ie)
211 area_nodes(7,nel_dg_loc) = alfa23(ie)
212
213 area_nodes(8,nel_dg_loc) = alfa31(ie)
214 area_nodes(9,nel_dg_loc) = alfa32(ie)
215 area_nodes(10,nel_dg_loc) = alfa33(ie)
216
217 area_nodes(11,nel_dg_loc) = beta11(ie)
218 area_nodes(12,nel_dg_loc) = beta12(ie)
219 area_nodes(13,nel_dg_loc) = beta13(ie)
220
221 area_nodes(14,nel_dg_loc) = beta21(ie)
222 area_nodes(15,nel_dg_loc) = beta22(ie)
223 area_nodes(16,nel_dg_loc) = beta23(ie)
224
225 area_nodes(17,nel_dg_loc) = beta31(ie)
226 area_nodes(18,nel_dg_loc) = beta32(ie)
227 area_nodes(19,nel_dg_loc) = beta33(ie)
228
229 area_nodes(20,nel_dg_loc) = gamma1(ie)
230 area_nodes(21,nel_dg_loc) = gamma2(ie)
231 area_nodes(22,nel_dg_loc) = gamma3(ie)
232
233 area_nodes(23,nel_dg_loc) = delta1(ie)
234 area_nodes(24,nel_dg_loc) = delta2(ie)
235 area_nodes(25,nel_dg_loc) = delta3(ie)
236
237 endif
238
239 !face z = - 1
240 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
241 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
242 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
243 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
244
245 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
246 nel_dg_loc = nel_dg_loc +1
247
248 call get_area_face(xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
249 zs(ne1), zs(ne2), zs(ne3), zs(ne4), surf, ie)
250
251 ! face z = 1
252
253 faces(1,nel_dg_loc) = tag_mat(im)
254 faces(2,nel_dg_loc) = local_el_num(ie)
255 faces(3,nel_dg_loc) = 3
256
257 area_nodes(1,nel_dg_loc) = surf
258 area_nodes(2,nel_dg_loc) = alfa11(ie)
259 area_nodes(3,nel_dg_loc) = alfa12(ie)
260 area_nodes(4,nel_dg_loc) = alfa13(ie)
261
262 area_nodes(5,nel_dg_loc) = alfa21(ie)
263 area_nodes(6,nel_dg_loc) = alfa22(ie)
264 area_nodes(7,nel_dg_loc) = alfa23(ie)
265
266 area_nodes(8,nel_dg_loc) = alfa31(ie)
267 area_nodes(9,nel_dg_loc) = alfa32(ie)
268 area_nodes(10,nel_dg_loc) = alfa33(ie)
269
270 area_nodes(11,nel_dg_loc) = beta11(ie)
271 area_nodes(12,nel_dg_loc) = beta12(ie)
272 area_nodes(13,nel_dg_loc) = beta13(ie)
273
274 area_nodes(14,nel_dg_loc) = beta21(ie)
275 area_nodes(15,nel_dg_loc) = beta22(ie)
276 area_nodes(16,nel_dg_loc) = beta23(ie)
277
278 area_nodes(17,nel_dg_loc) = beta31(ie)
279 area_nodes(18,nel_dg_loc) = beta32(ie)
280 area_nodes(19,nel_dg_loc) = beta33(ie)
281
282 area_nodes(20,nel_dg_loc) = gamma1(ie)
283 area_nodes(21,nel_dg_loc) = gamma2(ie)
284 area_nodes(22,nel_dg_loc) = gamma3(ie)
285
286 area_nodes(23,nel_dg_loc) = delta1(ie)
287 area_nodes(24,nel_dg_loc) = delta2(ie)
288 area_nodes(25,nel_dg_loc) = delta3(ie)
289
290
291 endif
292
293 ! face x = 1
294 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
295 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
296 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
297 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
298
299
300 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
301 nel_dg_loc = nel_dg_loc +1
302
303 call get_area_face(xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
304 zs(ne1), zs(ne2), zs(ne3), zs(ne4), surf, ie)
305
306 ! face x = - 1
307
308 faces(1,nel_dg_loc) = tag_mat(im)
309 faces(2,nel_dg_loc) = local_el_num(ie)
310 faces(3,nel_dg_loc) = 4
311
312 area_nodes(1,nel_dg_loc) = surf
313 area_nodes(2,nel_dg_loc) = alfa11(ie)
314 area_nodes(3,nel_dg_loc) = alfa12(ie)
315 area_nodes(4,nel_dg_loc) = alfa13(ie)
316
317 area_nodes(5,nel_dg_loc) = alfa21(ie)
318 area_nodes(6,nel_dg_loc) = alfa22(ie)
319 area_nodes(7,nel_dg_loc) = alfa23(ie)
320
321 area_nodes(8,nel_dg_loc) = alfa31(ie)
322 area_nodes(9,nel_dg_loc) = alfa32(ie)
323 area_nodes(10,nel_dg_loc) = alfa33(ie)
324
325 area_nodes(11,nel_dg_loc) = beta11(ie)
326 area_nodes(12,nel_dg_loc) = beta12(ie)
327 area_nodes(13,nel_dg_loc) = beta13(ie)
328
329 area_nodes(14,nel_dg_loc) = beta21(ie)
330 area_nodes(15,nel_dg_loc) = beta22(ie)
331 area_nodes(16,nel_dg_loc) = beta23(ie)
332
333 area_nodes(17,nel_dg_loc) = beta31(ie)
334 area_nodes(18,nel_dg_loc) = beta32(ie)
335 area_nodes(19,nel_dg_loc) = beta33(ie)
336
337 area_nodes(20,nel_dg_loc) = gamma1(ie)
338 area_nodes(21,nel_dg_loc) = gamma2(ie)
339 area_nodes(22,nel_dg_loc) = gamma3(ie)
340
341 area_nodes(23,nel_dg_loc) = delta1(ie)
342 area_nodes(24,nel_dg_loc) = delta2(ie)
343 area_nodes(25,nel_dg_loc) = delta3(ie)
344
345
346 endif
347
348 ! face y = 1
349 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
350 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
351 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
352 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
353
354
355 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
356 nel_dg_loc = nel_dg_loc +1
357
358 call get_area_face(xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
359 zs(ne1), zs(ne2), zs(ne3), zs(ne4), surf, ie)
360
361 faces(1,nel_dg_loc) = tag_mat(im)
362 faces(2,nel_dg_loc) = local_el_num(ie)
363 faces(3,nel_dg_loc) = 5
364
365 area_nodes(1,nel_dg_loc) = surf
366 area_nodes(2,nel_dg_loc) = alfa11(ie)
367 area_nodes(3,nel_dg_loc) = alfa12(ie)
368 area_nodes(4,nel_dg_loc) = alfa13(ie)
369
370 area_nodes(5,nel_dg_loc) = alfa21(ie)
371 area_nodes(6,nel_dg_loc) = alfa22(ie)
372 area_nodes(7,nel_dg_loc) = alfa23(ie)
373
374 area_nodes(8,nel_dg_loc) = alfa31(ie)
375 area_nodes(9,nel_dg_loc) = alfa32(ie)
376 area_nodes(10,nel_dg_loc) = alfa33(ie)
377
378 area_nodes(11,nel_dg_loc) = beta11(ie)
379 area_nodes(12,nel_dg_loc) = beta12(ie)
380 area_nodes(13,nel_dg_loc) = beta13(ie)
381
382 area_nodes(14,nel_dg_loc) = beta21(ie)
383 area_nodes(15,nel_dg_loc) = beta22(ie)
384 area_nodes(16,nel_dg_loc) = beta23(ie)
385
386 area_nodes(17,nel_dg_loc) = beta31(ie)
387 area_nodes(18,nel_dg_loc) = beta32(ie)
388 area_nodes(19,nel_dg_loc) = beta33(ie)
389
390 area_nodes(20,nel_dg_loc) = gamma1(ie)
391 area_nodes(21,nel_dg_loc) = gamma2(ie)
392 area_nodes(22,nel_dg_loc) = gamma3(ie)
393
394 area_nodes(23,nel_dg_loc) = delta1(ie)
395 area_nodes(24,nel_dg_loc) = delta2(ie)
396 area_nodes(25,nel_dg_loc) = delta3(ie)
397
398 endif
399
400 ! face z = 1
401 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
402 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
403 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
404 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
405
406
407 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
408 nel_dg_loc = nel_dg_loc +1
409
410 call get_area_face(xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
411 zs(ne1), zs(ne2), zs(ne3), zs(ne4), surf, ie)
412
413
414 !face z = - 1
415
416 faces(1,nel_dg_loc) = tag_mat(im)
417 faces(2,nel_dg_loc) = local_el_num(ie)
418 faces(3,nel_dg_loc) = 6
419
420 area_nodes(1,nel_dg_loc) = surf
421 area_nodes(2,nel_dg_loc) = alfa11(ie)
422 area_nodes(3,nel_dg_loc) = alfa12(ie)
423 area_nodes(4,nel_dg_loc) = alfa13(ie)
424
425 area_nodes(5,nel_dg_loc) = alfa21(ie)
426 area_nodes(6,nel_dg_loc) = alfa22(ie)
427 area_nodes(7,nel_dg_loc) = alfa23(ie)
428
429 area_nodes(8,nel_dg_loc) = alfa31(ie)
430 area_nodes(9,nel_dg_loc) = alfa32(ie)
431 area_nodes(10,nel_dg_loc) = alfa33(ie)
432
433 area_nodes(11,nel_dg_loc) = beta11(ie)
434 area_nodes(12,nel_dg_loc) = beta12(ie)
435 area_nodes(13,nel_dg_loc) = beta13(ie)
436
437 area_nodes(14,nel_dg_loc) = beta21(ie)
438 area_nodes(15,nel_dg_loc) = beta22(ie)
439 area_nodes(16,nel_dg_loc) = beta23(ie)
440
441 area_nodes(17,nel_dg_loc) = beta31(ie)
442 area_nodes(18,nel_dg_loc) = beta32(ie)
443 area_nodes(19,nel_dg_loc) = beta33(ie)
444
445 area_nodes(20,nel_dg_loc) = gamma1(ie)
446 area_nodes(21,nel_dg_loc) = gamma2(ie)
447 area_nodes(22,nel_dg_loc) = gamma3(ie)
448
449 area_nodes(23,nel_dg_loc) = delta1(ie)
450 area_nodes(24,nel_dg_loc) = delta2(ie)
451 area_nodes(25,nel_dg_loc) = delta3(ie)
452
453 endif
454 endif
455 enddo
456 enddo
457
458
459
460 return
461
462 end subroutine setup_dg
get_area_face
subroutine get_area_face(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, surf, ielem)
Computes area of a face (quad).
Definition GET_AREA_FACE.f90:39
setup_dg
subroutine setup_dg(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, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, faces, area_nodes)
Setup for DG faces.
Definition SETUP_DG.f90:115