SPEED
MAKE_DG_INTERFACE.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
77
78
79 subroutine make_dg_interface(nm, sd, tag_mat, prop_mat, cs_nnz_loc, cs_loc, &
80 nn_loc, local_n_num, ne_loc, local_el_num, &
81 xs, ys, zs,&
82 nel_dg_loc, nel_dg_glo, &
83 i4count, mpi_id, mpi_comm, mpi_np,&
84 alfa11, alfa12, alfa13, &
85 alfa21, alfa22, alfa23, &
86 alfa31, alfa32, alfa33, &
87 beta11, beta12, beta13, &
88 beta21, beta22, beta23, &
89 beta31, beta32, beta33, &
90 gamma1, gamma2, gamma3, &
91 delta1, delta2, delta3, dg_cnst, penalty_c, &
92 faces, area_nodes, el_new, filename, testmode)
93
94
95 use max_var
96 use str_mesh_after
97 use dgjump
98
99 implicit none
100
101 type(element_after), dimension(:), allocatable :: dg_els
102 type(el4loop), dimension(nel_dg_loc), intent(inout):: el_new
103
104 character*70 :: filename, filename1, filename2
105
106 integer*4 :: mpi_comm, mpi_id, mpi_ierr, mpi_np, ishift, jshift, nel_frc
107 integer*4 :: im, nn, ie, ned, ip, mm, nnz, nnz_p, nnz_m, nnz_p_only_uv, nnz_m_only_uv
108 integer*4 :: ne1, ne2, ne3, ne4, ic1, ic2, ic3, ic4, ifr, fracture_yn
109 integer*4 :: ne5, ne6, ne7, ne8, ic5, ic6, ic7, ic8
110 integer*4 :: nm, cs_nnz_loc, nn_loc, ne_loc, nel_dg_loc, nel_dg_glo
111 integer*4 :: n_line, ielem, iface, iene, ifacene, imne
112 integer*4 :: int_trash, statuss, i, tt, ih, it, p, j,k
113 integer*4 :: unit_interface, unitname, unitname1, unitname2, unitname_frc
114 integer*4 :: error
115 integer*4 :: mpierror
116 integer*4 :: nofne_el, ic, testmode
117
118 integer*4, dimension(nm) :: tag_mat, sd
119 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
120 integer*4, dimension(nn_loc) :: local_n_num, i4count
121 integer*4, dimension(ne_loc) :: local_el_num
122 integer*4, dimension(:), allocatable :: I4S, J4S
123
124 integer*4, dimension(:,:), allocatable :: con_DG, con_DG_frc_int
125 integer*4, dimension(:,:), allocatable :: copia
126 integer*4, dimension(3,nel_dg_glo) :: faces
127
128 real*8 :: real_trash, lambda, mu, pen_p, pen_h, pen, penalty_c, dg_cnst, dummy
129 real*8 :: cp_a, cp_b, cp_c, cp_d, cp_e, cp_f, cp_g, cp_n, cp_p
130 real*8 :: csi, eta, zeta, normal_x, normal_y, normal_z
131 real*8 :: c_alfa11, c_alfa12, c_alfa13, c_alfa21, c_alfa22, c_alfa23, c_alfa31, c_alfa32, c_alfa33
132 real*8 :: c_beta11, c_beta12, c_beta13, c_beta21, c_beta22, c_beta23, c_beta31, c_beta32, c_beta33
133 real*8 :: c_gamma1, c_gamma2, c_gamma3, c_delta1, c_delta2, c_delta3
134 real*8 :: z11,z12,z13,z22,z21,z23,z31,z32,z33, detz, zn_coef, zt_coef
135
136
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
146 real*8, dimension(:), allocatable :: m4s
147
148 real*8, dimension(:,:), allocatable :: nodes_dg, con_dg_frc_real
149 real*8, dimension(:,:), allocatable :: jp,jm, jp_only_uv,jm_only_uv
150 real*8, dimension(nn_loc) :: xs,ys,zs
151 real*8, dimension(nm,4) :: prop_mat
152 real*8, dimension(25,nel_dg_glo) :: area_nodes
153 real*8, dimension(3,3) :: invz, compz, ide
154
155
156
157
158! Reading DGFS.input file'
159 unitname = 50
160 open(unitname,file=filename)
161 read(unitname,*) n_line
162
163 allocate(con_dg(6,2*n_line), nodes_dg(6,2*n_line))
164 con_dg = 0
165 nodes_dg = 0.d0
166
167
168 do i = 1, n_line
169 read(unitname,*) con_dg(1,i), con_dg(2,i), con_dg(3,i), con_dg(4,i), con_dg(5,i), con_dg(6,i), &
170 nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), nodes_dg(4,i), nodes_dg(5,i), nodes_dg(6,i)
171 enddo
172
173
174 close(unitname)
175
176
177 do i = 1, 3
178 con_dg(i, n_line+1: 2*n_line) = con_dg(3+i,1:n_line)
179 enddo
180
181 do i = 1, 3
182 con_dg(3+i, n_line+1: 2*n_line) = con_dg(i,1:n_line)
183 enddo
184
185
186
187 do i = 1, 6
188 nodes_dg(i, n_line+1: 2*n_line) = nodes_dg(i,1:n_line)
189 enddo
190
191
192! Reading NORMALL.input file for fracture
193
194 unitname = 50
195 open(unitname,file='NORMALL.input')
196 read(unitname,*) nel_frc
197
198 allocate(con_dg_frc_int(4,nel_frc), con_dg_frc_real(2,nel_frc))
199 con_dg_frc_real = 0; con_dg_frc_int = 0;
200
201 do i = 1, nel_frc
202 read(unitname,*) con_dg_frc_int(1,i), con_dg_frc_int(2,i), con_dg_frc_int(3,i),&
203 dummy, dummy, dummy, con_dg_frc_int(4,i), con_dg_frc_real(1,i), con_dg_frc_real(2,i)
204 enddo
205
206
207 close(unitname)
208
209
210
211 allocate(dg_els(nel_dg_loc))
212
213!***************************************************************************************************************
214
215 nel_dg_loc = 0
216 ned = cs_loc(0) - 1
217
218
219 do im = 1,nm
220
221 nn = sd(im) +1
222
223 do ie = 1,ned
224 if (cs_loc(cs_loc(ie -1) + 0) .eq. tag_mat(im)) then
225
226 ! face 1: x = - 1
227 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
228 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
229 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
230 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
231
232
233 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
234
235 ip = 0
236 nel_dg_loc = nel_dg_loc + 1
237 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
238 dg_els(nel_dg_loc)%face_el = 1
239 dg_els(nel_dg_loc)%mat = tag_mat(im)
240 dg_els(nel_dg_loc)%spct_deg = nn-1
241 dg_els(nel_dg_loc)%quad_rule = ip
242
243 i = 1
244 do while (i .le. 2*n_line)
245
246 if(con_dg(2,i) .eq. local_el_num(ie) .and. con_dg(3,i) .eq. 1) then
247
248
249 ip = ip + 1
250 dg_els(nel_dg_loc)%quad_rule = ip
251
252 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
253 alfa11(ie),alfa12(ie),alfa13(ie), &
254 alfa21(ie),alfa22(ie),alfa23(ie), &
255 alfa31(ie),alfa32(ie),alfa33(ie), &
256 beta11(ie),beta12(ie),beta13(ie), &
257 beta21(ie),beta22(ie),beta23(ie), &
258 beta31(ie),beta32(ie),beta33(ie), &
259 gamma1(ie),gamma2(ie),gamma3(ie), &
260 delta1(ie),delta2(ie),delta3(ie), &
261 tt, csi, eta, zeta, nofinr, mpi_id,&
262 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
263
264
265 dg_els(nel_dg_loc)%x_pl(ip) = -1.d0
266 dg_els(nel_dg_loc)%y_pl(ip) = eta
267 dg_els(nel_dg_loc)%z_pl(ip) = zeta
268 dg_els(nel_dg_loc)%wx_pl(ip) = nodes_dg(4,i)
269 dg_els(nel_dg_loc)%wy_pl(ip) = nodes_dg(5,i)
270 dg_els(nel_dg_loc)%wz_pl(ip) = nodes_dg(6,i)
271
272 dg_els(nel_dg_loc)%omega_minus(ip,0) = con_dg(4,i)
273 dg_els(nel_dg_loc)%omega_minus(ip,1) = con_dg(5,i)
274 dg_els(nel_dg_loc)%omega_minus(ip,2) = con_dg(6,i)
275 dg_els(nel_dg_loc)%omega_minus(ip,3) = 0
276
277
278 ! find the corresponding neigh. elem in faces matrix
279
280 call get_face_dg(faces, nel_dg_glo, con_dg(5,i), con_dg(6,i), ih)
281
282 call make_bilinear_map(area_nodes(2:25,ih), &
283 c_alfa11, c_alfa12, c_alfa13, &
284 c_alfa21, c_alfa22, c_alfa23, &
285 c_alfa31, c_alfa32, c_alfa33, &
286 c_beta11, c_beta12, c_beta13, &
287 c_beta21, c_beta22, c_beta23, &
288 c_beta31, c_beta32, c_beta33, &
289 c_gamma1, c_gamma2, c_gamma3, &
290 c_delta1, c_delta2, c_delta3)
291
292
293 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
294 c_alfa11, c_alfa12, c_alfa13, &
295 c_alfa21, c_alfa22, c_alfa23, &
296 c_alfa31, c_alfa32, c_alfa33, &
297 c_beta11, c_beta12, c_beta13, &
298 c_beta21, c_beta22, c_beta23, &
299 c_beta31, c_beta32, c_beta33, &
300 c_gamma1, c_gamma2, c_gamma3, &
301 c_delta1, c_delta2, c_delta3, &
302 tt, csi, eta, zeta, nofinr, mpi_id,&
303 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
304
305 select case(con_dg(6,i))
306 case(1)
307 csi = -1.d0
308 case(2)
309 eta = -1.d0
310 case(3)
311 zeta = -1.d0
312 case(4)
313 csi = 1.d0
314 case(5)
315 eta = 1.d0
316 case(6)
317 zeta = 1.d0
318 case default
319 write(*,*) 'error in make_interface'
320 end select
321
322 dg_els(nel_dg_loc)%x_mn(ip) = csi
323 dg_els(nel_dg_loc)%y_mn(ip) = eta
324 dg_els(nel_dg_loc)%z_mn(ip) = zeta
325
326
327 endif
328 i = i + 1
329
330 enddo
331
332
333
334
335 call make_normal(1,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
336 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
337
338 dg_els(nel_dg_loc)%nx = normal_x
339 dg_els(nel_dg_loc)%ny = normal_y
340 dg_els(nel_dg_loc)%nz = normal_z
341
342
343
344
345 endif
346
347
348 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
349 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
350 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
351 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
352
353
354 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
355
356 !face 2: y = -1
357 ip = 0
358 nel_dg_loc = nel_dg_loc +1
359 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
360 dg_els(nel_dg_loc)%face_el = 2
361 dg_els(nel_dg_loc)%mat = tag_mat(im)
362 dg_els(nel_dg_loc)%spct_deg = nn-1
363 dg_els(nel_dg_loc)%quad_rule = ip
364 i = 1
365
366
367 do while (i .le. 2*n_line)
368
369 if(con_dg(2,i) .eq. local_el_num(ie) .and. con_dg(3,i) .eq. 2) then
370
371 ip = ip + 1
372 dg_els(nel_dg_loc)%quad_rule = ip
373
374 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
375 alfa11(ie),alfa12(ie),alfa13(ie), &
376 alfa21(ie),alfa22(ie),alfa23(ie), &
377 alfa31(ie),alfa32(ie),alfa33(ie), &
378 beta11(ie),beta12(ie),beta13(ie), &
379 beta21(ie),beta22(ie),beta23(ie), &
380 beta31(ie),beta32(ie),beta33(ie), &
381 gamma1(ie),gamma2(ie),gamma3(ie), &
382 delta1(ie),delta2(ie),delta3(ie), &
383 tt, csi, eta, zeta, nofinr, mpi_id,&
384 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
385
386
387
388 dg_els(nel_dg_loc)%x_pl(ip) = csi
389 dg_els(nel_dg_loc)%y_pl(ip) = -1.d0
390 dg_els(nel_dg_loc)%z_pl(ip) = zeta
391 dg_els(nel_dg_loc)%wx_pl(ip) = nodes_dg(4,i)
392 dg_els(nel_dg_loc)%wy_pl(ip) = nodes_dg(5,i)
393 dg_els(nel_dg_loc)%wz_pl(ip) = nodes_dg(6,i)
394
395 dg_els(nel_dg_loc)%omega_minus(ip,0) = con_dg(4,i)
396 dg_els(nel_dg_loc)%omega_minus(ip,1) = con_dg(5,i)
397 dg_els(nel_dg_loc)%omega_minus(ip,2) = con_dg(6,i)
398 dg_els(nel_dg_loc)%omega_minus(ip,3) = 0
399
400
401 ! find the corresponding neigh. elem in faces matrix
402
403 call get_face_dg(faces, nel_dg_glo, con_dg(5,i), con_dg(6,i), ih)
404
405 call make_bilinear_map(area_nodes(2:25,ih), &
406 c_alfa11, c_alfa12, c_alfa13, &
407 c_alfa21, c_alfa22, c_alfa23, &
408 c_alfa31, c_alfa32, c_alfa33, &
409 c_beta11, c_beta12, c_beta13, &
410 c_beta21, c_beta22, c_beta23, &
411 c_beta31, c_beta32, c_beta33, &
412 c_gamma1, c_gamma2, c_gamma3, &
413 c_delta1, c_delta2, c_delta3)
414
415
416 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
417 c_alfa11, c_alfa12, c_alfa13, &
418 c_alfa21, c_alfa22, c_alfa23, &
419 c_alfa31, c_alfa32, c_alfa33, &
420 c_beta11, c_beta12, c_beta13, &
421 c_beta21, c_beta22, c_beta23, &
422 c_beta31, c_beta32, c_beta33, &
423 c_gamma1, c_gamma2, c_gamma3, &
424 c_delta1, c_delta2, c_delta3, &
425 tt, csi, eta, zeta, nofinr, mpi_id,&
426 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
427
428 select case(con_dg(6,i))
429 case(1)
430 csi = -1.d0
431 case(2)
432 eta = -1.d0
433 case(3)
434 zeta = -1.d0
435 case(4)
436 csi = 1.d0
437 case(5)
438 eta = 1.d0
439 case(6)
440 zeta = 1.d0
441 case default
442 write(*,*) 'error in make_interface'
443 end select
444
445 dg_els(nel_dg_loc)%x_mn(ip) = csi
446 dg_els(nel_dg_loc)%y_mn(ip) = eta
447 dg_els(nel_dg_loc)%z_mn(ip) = zeta
448
449 endif
450
451 i = i+1
452
453 enddo
454
455 call make_normal(2,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
456 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
457
458 dg_els(nel_dg_loc)%nx = normal_x
459 dg_els(nel_dg_loc)%ny = normal_y
460 dg_els(nel_dg_loc)%nz = normal_z
461
462 endif
463
464
465 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
466 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
467 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
468 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
469
470 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
471
472 !face 3: z = -1
473 ip = 0
474 nel_dg_loc = nel_dg_loc +1
475 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
476 dg_els(nel_dg_loc)%face_el = 3
477 dg_els(nel_dg_loc)%mat = tag_mat(im)
478 dg_els(nel_dg_loc)%spct_deg = nn-1
479 dg_els(nel_dg_loc)%quad_rule = ip
480
481 i = 1
482
483 do while (i .le. 2*n_line)
484
485 if(con_dg(2,i) .eq. local_el_num(ie) .and. con_dg(3,i) .eq. 3) then
486
487
488 ip = ip + 1
489 dg_els(nel_dg_loc)%quad_rule = ip
490
491 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
492 alfa11(ie),alfa12(ie),alfa13(ie), &
493 alfa21(ie),alfa22(ie),alfa23(ie), &
494 alfa31(ie),alfa32(ie),alfa33(ie), &
495 beta11(ie),beta12(ie),beta13(ie), &
496 beta21(ie),beta22(ie),beta23(ie), &
497 beta31(ie),beta32(ie),beta33(ie), &
498 gamma1(ie),gamma2(ie),gamma3(ie), &
499 delta1(ie),delta2(ie),delta3(ie), &
500 tt, csi, eta, zeta, nofinr, mpi_id,&
501 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
502
503
504
505 dg_els(nel_dg_loc)%x_pl(ip) = csi
506 dg_els(nel_dg_loc)%y_pl(ip) = eta
507 dg_els(nel_dg_loc)%z_pl(ip) = -1.d0
508 dg_els(nel_dg_loc)%wx_pl(ip) = nodes_dg(4,i)
509 dg_els(nel_dg_loc)%wy_pl(ip) = nodes_dg(5,i)
510 dg_els(nel_dg_loc)%wz_pl(ip) = nodes_dg(6,i)
511
512 dg_els(nel_dg_loc)%omega_minus(ip,0) = con_dg(4,i)
513 dg_els(nel_dg_loc)%omega_minus(ip,1) = con_dg(5,i)
514 dg_els(nel_dg_loc)%omega_minus(ip,2) = con_dg(6,i)
515 dg_els(nel_dg_loc)%omega_minus(ip,3) = 0
516
517
518 ! find the corresponding neigh. elem in faces matrix
519
520 call get_face_dg(faces, nel_dg_glo, con_dg(5,i), con_dg(6,i), ih)
521
522 call make_bilinear_map(area_nodes(2:25,ih), &
523 c_alfa11, c_alfa12, c_alfa13, &
524 c_alfa21, c_alfa22, c_alfa23, &
525 c_alfa31, c_alfa32, c_alfa33, &
526 c_beta11, c_beta12, c_beta13, &
527 c_beta21, c_beta22, c_beta23, &
528 c_beta31, c_beta32, c_beta33, &
529 c_gamma1, c_gamma2, c_gamma3, &
530 c_delta1, c_delta2, c_delta3)
531
532
533 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
534 c_alfa11, c_alfa12, c_alfa13, &
535 c_alfa21, c_alfa22, c_alfa23, &
536 c_alfa31, c_alfa32, c_alfa33, &
537 c_beta11, c_beta12, c_beta13, &
538 c_beta21, c_beta22, c_beta23, &
539 c_beta31, c_beta32, c_beta33, &
540 c_gamma1, c_gamma2, c_gamma3, &
541 c_delta1, c_delta2, c_delta3, &
542 tt, csi, eta, zeta, nofinr, mpi_id,&
543 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
544
545 select case(con_dg(6,i))
546 case(1)
547 csi = -1.d0
548 case(2)
549 eta = -1.d0
550 case(3)
551 zeta = -1.d0
552 case(4)
553 csi = 1.d0
554 case(5)
555 eta = 1.d0
556 case(6)
557 zeta = 1.d0
558 case default
559 write(*,*) 'error in make_interface'
560 end select
561
562 dg_els(nel_dg_loc)%x_mn(ip) = csi
563 dg_els(nel_dg_loc)%y_mn(ip) = eta
564 dg_els(nel_dg_loc)%z_mn(ip) = zeta
565
566 endif
567 i = i +1
568
569
570 enddo
571
572 call make_normal(3,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
573 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
574
575 dg_els(nel_dg_loc)%nx = normal_x
576 dg_els(nel_dg_loc)%ny = normal_y
577 dg_els(nel_dg_loc)%nz = normal_z
578
579 endif
580
581 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
582 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
583 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
584 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
585
586
587 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
588
589 !face 4 : x = 1
590 ip = 0
591 nel_dg_loc = nel_dg_loc +1
592 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
593 dg_els(nel_dg_loc)%face_el = 4
594 dg_els(nel_dg_loc)%mat = tag_mat(im)
595 dg_els(nel_dg_loc)%spct_deg = nn-1
596 dg_els(nel_dg_loc)%quad_rule = ip
597 i = 1
598
599 do while (i .le. 2*n_line)
600
601 if(con_dg(2,i) .eq. local_el_num(ie) .and. con_dg(3,i) .eq. 4) then
602
603
604 ip = ip + 1
605 dg_els(nel_dg_loc)%quad_rule = ip
606
607 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
608 alfa11(ie),alfa12(ie),alfa13(ie), &
609 alfa21(ie),alfa22(ie),alfa23(ie), &
610 alfa31(ie),alfa32(ie),alfa33(ie), &
611 beta11(ie),beta12(ie),beta13(ie), &
612 beta21(ie),beta22(ie),beta23(ie), &
613 beta31(ie),beta32(ie),beta33(ie), &
614 gamma1(ie),gamma2(ie),gamma3(ie), &
615 delta1(ie),delta2(ie),delta3(ie), &
616 tt, csi, eta, zeta, nofinr, mpi_id,&
617 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
618
619 dg_els(nel_dg_loc)%x_pl(ip) = 1.d0
620 dg_els(nel_dg_loc)%y_pl(ip) = eta
621 dg_els(nel_dg_loc)%z_pl(ip) = zeta
622 dg_els(nel_dg_loc)%wx_pl(ip) = nodes_dg(4,i)
623 dg_els(nel_dg_loc)%wy_pl(ip) = nodes_dg(5,i)
624 dg_els(nel_dg_loc)%wz_pl(ip) = nodes_dg(6,i)
625
626 dg_els(nel_dg_loc)%omega_minus(ip,0) = con_dg(4,i)
627 dg_els(nel_dg_loc)%omega_minus(ip,1) = con_dg(5,i)
628 dg_els(nel_dg_loc)%omega_minus(ip,2) = con_dg(6,i)
629 dg_els(nel_dg_loc)%omega_minus(ip,3) = 0
630
631
632 ! find the corresponding neigh. elem in faces matrix
633
634 call get_face_dg(faces, nel_dg_glo, con_dg(5,i), con_dg(6,i), ih)
635
636 call make_bilinear_map(area_nodes(2:25,ih), &
637 c_alfa11, c_alfa12, c_alfa13, &
638 c_alfa21, c_alfa22, c_alfa23, &
639 c_alfa31, c_alfa32, c_alfa33, &
640 c_beta11, c_beta12, c_beta13, &
641 c_beta21, c_beta22, c_beta23, &
642 c_beta31, c_beta32, c_beta33, &
643 c_gamma1, c_gamma2, c_gamma3, &
644 c_delta1, c_delta2, c_delta3)
645
646
647 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
648 c_alfa11, c_alfa12, c_alfa13, &
649 c_alfa21, c_alfa22, c_alfa23, &
650 c_alfa31, c_alfa32, c_alfa33, &
651 c_beta11, c_beta12, c_beta13, &
652 c_beta21, c_beta22, c_beta23, &
653 c_beta31, c_beta32, c_beta33, &
654 c_gamma1, c_gamma2, c_gamma3, &
655 c_delta1, c_delta2, c_delta3, &
656 tt, csi, eta, zeta, nofinr, mpi_id,&
657 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
658
659 select case(con_dg(6,i))
660 case(1)
661 csi = -1.d0
662 case(2)
663 eta = -1.d0
664 case(3)
665 zeta = -1.d0
666 case(4)
667 csi = 1.d0
668 case(5)
669 eta = 1.d0
670 case(6)
671 zeta = 1.d0
672 case default
673 write(*,*) 'error in make_interface'
674 end select
675
676 dg_els(nel_dg_loc)%x_mn(ip) = csi
677 dg_els(nel_dg_loc)%y_mn(ip) = eta
678 dg_els(nel_dg_loc)%z_mn(ip) = zeta
679
680 endif
681 i = i+1
682
683 enddo
684
685 call make_normal(4,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
686 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
687
688 dg_els(nel_dg_loc)%nx = normal_x
689 dg_els(nel_dg_loc)%ny = normal_y
690 dg_els(nel_dg_loc)%nz = normal_z
691
692 endif
693
694
695 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
696 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
697 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
698 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
699
700
701 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
702
703 !face 5: y = 1
704 ip = 0
705 nel_dg_loc = nel_dg_loc +1
706 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
707 dg_els(nel_dg_loc)%face_el = 5
708 dg_els(nel_dg_loc)%mat = tag_mat(im)
709 dg_els(nel_dg_loc)%spct_deg = nn-1
710 dg_els(nel_dg_loc)%quad_rule = ip
711 i = 1
712
713 do while (i .le. 2*n_line)
714
715
716 if(con_dg(2,i) .eq. local_el_num(ie) .and. con_dg(3,i) .eq. 5) then
717
718 ip = ip + 1
719 dg_els(nel_dg_loc)%quad_rule = ip
720
721 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
722 alfa11(ie),alfa12(ie),alfa13(ie), &
723 alfa21(ie),alfa22(ie),alfa23(ie), &
724 alfa31(ie),alfa32(ie),alfa33(ie), &
725 beta11(ie),beta12(ie),beta13(ie), &
726 beta21(ie),beta22(ie),beta23(ie), &
727 beta31(ie),beta32(ie),beta33(ie), &
728 gamma1(ie),gamma2(ie),gamma3(ie), &
729 delta1(ie),delta2(ie),delta3(ie), &
730 tt, csi, eta, zeta, nofinr, mpi_id,&
731 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
732
733
734
735 dg_els(nel_dg_loc)%x_pl(ip) = csi
736 dg_els(nel_dg_loc)%y_pl(ip) = 1.d0
737 dg_els(nel_dg_loc)%z_pl(ip) = zeta
738 dg_els(nel_dg_loc)%wx_pl(ip) = nodes_dg(4,i)
739 dg_els(nel_dg_loc)%wy_pl(ip) = nodes_dg(5,i)
740 dg_els(nel_dg_loc)%wz_pl(ip) = nodes_dg(6,i)
741
742 dg_els(nel_dg_loc)%omega_minus(ip,0) = con_dg(4,i)
743 dg_els(nel_dg_loc)%omega_minus(ip,1) = con_dg(5,i)
744 dg_els(nel_dg_loc)%omega_minus(ip,2) = con_dg(6,i)
745 dg_els(nel_dg_loc)%omega_minus(ip,3) = 0
746
747
748 ! find the corresponding neigh. elem in faces matrix
749
750 call get_face_dg(faces, nel_dg_glo, con_dg(5,i), con_dg(6,i), ih)
751
752 call make_bilinear_map(area_nodes(2:25,ih), &
753 c_alfa11, c_alfa12, c_alfa13, &
754 c_alfa21, c_alfa22, c_alfa23, &
755 c_alfa31, c_alfa32, c_alfa33, &
756 c_beta11, c_beta12, c_beta13, &
757 c_beta21, c_beta22, c_beta23, &
758 c_beta31, c_beta32, c_beta33, &
759 c_gamma1, c_gamma2, c_gamma3, &
760 c_delta1, c_delta2, c_delta3)
761
762
763 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
764 c_alfa11, c_alfa12, c_alfa13, &
765 c_alfa21, c_alfa22, c_alfa23, &
766 c_alfa31, c_alfa32, c_alfa33, &
767 c_beta11, c_beta12, c_beta13, &
768 c_beta21, c_beta22, c_beta23, &
769 c_beta31, c_beta32, c_beta33, &
770 c_gamma1, c_gamma2, c_gamma3, &
771 c_delta1, c_delta2, c_delta3, &
772 tt, csi, eta, zeta, nofinr, mpi_id,&
773 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
774
775 select case(con_dg(6,i))
776 case(1)
777 csi = -1.d0
778 case(2)
779 eta = -1.d0
780 case(3)
781 zeta = -1.d0
782 case(4)
783 csi = 1.d0
784 case(5)
785 eta = 1.d0
786 case(6)
787 zeta = 1.d0
788 case default
789 write(*,*) 'error in make_interface'
790 end select
791
792 dg_els(nel_dg_loc)%x_mn(ip) = csi
793 dg_els(nel_dg_loc)%y_mn(ip) = eta
794 dg_els(nel_dg_loc)%z_mn(ip) = zeta
795
796 endif
797
798 i = i+1
799
800
801 enddo
802
803 call make_normal(5,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
804 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
805
806 dg_els(nel_dg_loc)%nx = normal_x
807 dg_els(nel_dg_loc)%ny = normal_y
808 dg_els(nel_dg_loc)%nz = normal_z
809
810
811 endif
812
813 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
814 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
815 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
816 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
817
818
819 if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
820
821 !face 6: z = 1
822 ip = 0
823 nel_dg_loc = nel_dg_loc + 1
824 dg_els(nel_dg_loc)%ind_el = local_el_num(ie)
825 dg_els(nel_dg_loc)%face_el = 6
826 dg_els(nel_dg_loc)%mat = tag_mat(im)
827 dg_els(nel_dg_loc)%spct_deg = nn-1
828 dg_els(nel_dg_loc)%quad_rule = ip
829 i = 1
830
831 do while (i .le. 2*n_line)
832
833 if(con_dg(2,i) .eq. local_el_num(ie) .and. con_dg(3,i) .eq. 6) then
834
835 ip = ip + 1
836 dg_els(nel_dg_loc)%quad_rule = ip
837
838 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
839 alfa11(ie),alfa12(ie),alfa13(ie), &
840 alfa21(ie),alfa22(ie),alfa23(ie), &
841 alfa31(ie),alfa32(ie),alfa33(ie), &
842 beta11(ie),beta12(ie),beta13(ie), &
843 beta21(ie),beta22(ie),beta23(ie), &
844 beta31(ie),beta32(ie),beta33(ie), &
845 gamma1(ie),gamma2(ie),gamma3(ie), &
846 delta1(ie),delta2(ie),delta3(ie), &
847 tt, csi, eta, zeta, nofinr, mpi_id,&
848 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
849
850
851
852 dg_els(nel_dg_loc)%x_pl(ip) = csi
853 dg_els(nel_dg_loc)%y_pl(ip) = eta
854 dg_els(nel_dg_loc)%z_pl(ip) = 1.d0
855 dg_els(nel_dg_loc)%wx_pl(ip) = nodes_dg(4,i)
856 dg_els(nel_dg_loc)%wy_pl(ip) = nodes_dg(5,i)
857 dg_els(nel_dg_loc)%wz_pl(ip) = nodes_dg(6,i)
858
859
860
861 dg_els(nel_dg_loc)%omega_minus(ip,0) = con_dg(4,i)
862 dg_els(nel_dg_loc)%omega_minus(ip,1) = con_dg(5,i)
863 dg_els(nel_dg_loc)%omega_minus(ip,2) = con_dg(6,i)
864 dg_els(nel_dg_loc)%omega_minus(ip,3) = 0
865
866
867 ! find the corresponding neigh. elem in faces matrix
868
869 call get_face_dg(faces, nel_dg_glo, con_dg(5,i), con_dg(6,i), ih)
870
871 call make_bilinear_map(area_nodes(2:25,ih), &
872 c_alfa11, c_alfa12, c_alfa13, &
873 c_alfa21, c_alfa22, c_alfa23, &
874 c_alfa31, c_alfa32, c_alfa33, &
875 c_beta11, c_beta12, c_beta13, &
876 c_beta21, c_beta22, c_beta23, &
877 c_beta31, c_beta32, c_beta33, &
878 c_gamma1, c_gamma2, c_gamma3, &
879 c_delta1, c_delta2, c_delta3)
880
881
882 call newton_rapson(nodes_dg(1,i), nodes_dg(2,i), nodes_dg(3,i), &
883 c_alfa11, c_alfa12, c_alfa13, &
884 c_alfa21, c_alfa22, c_alfa23, &
885 c_alfa31, c_alfa32, c_alfa33, &
886 c_beta11, c_beta12, c_beta13, &
887 c_beta21, c_beta22, c_beta23, &
888 c_beta31, c_beta32, c_beta33, &
889 c_gamma1, c_gamma2, c_gamma3, &
890 c_delta1, c_delta2, c_delta3, &
891 tt, csi, eta, zeta, nofinr, mpi_id,&
892 con_dg(2,i),con_dg(5,i), 1.d-6, 1.01d0,1)
893
894 select case(con_dg(6,i))
895 case(1)
896 csi = -1.d0
897 case(2)
898 eta = -1.d0
899 case(3)
900 zeta = -1.d0
901 case(4)
902 csi = 1.d0
903 case(5)
904 eta = 1.d0
905 case(6)
906 zeta = 1.d0
907 case default
908 write(*,*) 'error in make_interface'
909 end select
910
911 dg_els(nel_dg_loc)%x_mn(ip) = csi
912 dg_els(nel_dg_loc)%y_mn(ip) = eta
913 dg_els(nel_dg_loc)%z_mn(ip) = zeta
914
915 endif
916
917 i = i+1
918 enddo
919
920 call make_normal(6,xs(ne1), xs(ne2), xs(ne3), xs(ne4), ys(ne1), ys(ne2), ys(ne3), ys(ne4), &
921 zs(ne1), zs(ne2), zs(ne3), zs(ne4), normal_x, normal_y, normal_z, 0, 0)
922
923 dg_els(nel_dg_loc)%nx = normal_x
924 dg_els(nel_dg_loc)%ny = normal_y
925 dg_els(nel_dg_loc)%nz = normal_z
926
927
928 endif
929 endif
930 enddo
931 enddo
932
933
934 deallocate(con_dg, nodes_dg)
935
936
937! Finding/ordering neigh. elements
938
939 do it = 1, nel_dg_loc
940 allocate(copia(max_quad_points,0:3))
941 copia = dg_els(it)%omega_minus
942
943 call get_neighbour_elem(copia, dg_els(it)%quad_rule, nofne_el, mpi_id, max_quad_points)
944
945 dg_els(it)%nofne = nofne_el
946 dg_els(it)%omega_minus(:,3) = copia(:,3)
947
948 deallocate(copia)
949 enddo
950
951
952!CHECK IF ALL DG EL HAVE A NEIGHBOUR
953 do it = 1, nel_dg_loc
954 if(dg_els(it)%nofne == 0) then
955 write(*,*) 'ATTENTION : element ', dg_els(it)%ind_el, ' has 0 neighbouring elements'
956 endif
957 enddo
958
959
960 do it = 1, nel_dg_loc
961
962 ic = 1
963 do while(ic .le. dg_els(it)%nofne)
964 do ip = 1, (dg_els(it)%quad_rule)
965 if(dg_els(it)%omega_minus(ip,3) .eq. ic) then
966 dg_els(it)%conf(ic,0) = dg_els(it)%omega_minus(ip,0)
967 dg_els(it)%conf(ic,1) = dg_els(it)%omega_minus(ip,1)
968 dg_els(it)%conf(ic,2) = dg_els(it)%omega_minus(ip,2)
969
970 endif
971 enddo
972 ic = ic + 1
973 enddo
974 enddo
975
976
977!*************************************************************************************************************
978! JUMP FOR DG INTERFACES
979!*************************************************************************************************************
980
981 !Needs OpenMP for speed-up
982 !Remove OpenMP from MAKE_LOC_MATRIX_DG.f90
983
984!!$OMP PARALLEL &
985!!$OMP PRIVATE(ie, nnz_p, nnz_m, ielem, iface, nn, im, ic1, ne1, ic) &
986!!$OMP PRIVATE(imne, iene, ifacene, mm, i, lambda, mu, ic2, JP, JM, pen_h, pen_p, pen) &
987!!$OMP PRIVATE(cp_a, cp_b, cp_c, cp_d, cp_e, cp_f, cp_g, cp_n, cp_p, mpi_id, j) &
988!!$OMP PRIVATE(c_alfa11, c_alfa12, c_alfa13, c_alfa21, c_alfa22, c_alfa23) &
989!!$OMP PRIVATE(c_alfa31, c_alfa32, c_alfa33, c_beta11, c_beta12, c_beta13) &
990!!$OMP PRIVATE(c_beta21, c_beta22, c_beta23, c_beta31, c_beta32, c_beta33) &
991!!$OMP PRIVATE(c_gamma1, c_gamma2, c_gamma3, c_delta1, c_delta2, c_delta3) &
992!!$OMP PRIVATE(k, I4S, ishift, jshift, J4S, M4S)
993
994!!$OMP DO
995
996 do ie = 1, nel_dg_loc
997 nnz_p = 0
998 nnz_m = 0
999 el_new(ie)%nnz_col = 0
1000 if(testmode .eq. 1) then
1001 el_new(ie)%nnz_col_only_uv = 0; nnz_p_only_uv = 0; nnz_m_only_uv = 0
1002 endif
1003
1004 ielem = dg_els(ie)%ind_el
1005 iface = dg_els(ie)%face_el
1006 nn = dg_els(ie)%spct_deg + 1
1007 im = dg_els(ie)%mat
1008
1009 el_new(ie)%ind = ielem
1010 el_new(ie)%face = iface
1011 el_new(ie)%deg = nn
1012 el_new(ie)%mate = im
1013
1014 el_new(ie)%num_of_ne = dg_els(ie)%nofne
1015
1016 call get_face_dg(faces, nel_dg_glo, ielem, iface, ic1)
1017 call get_indloc_from_indglo(local_el_num, ne_loc, ielem, ne1)
1018
1019 call get_face_dg(con_dg_frc_int(1:3,:), nel_dg_glo, ielem, iface, ifr)
1020 fracture_yn = con_dg_frc_int(4,ifr)
1021 zt_coef = con_dg_frc_real(1,ifr)
1022 zn_coef = con_dg_frc_real(2,ifr)
1023 compz = 0.d0; invz = 0.d0;
1024
1025 if(fracture_yn .eq. 1) then
1026 compz(1,1) = zn_coef*dg_els(ie)%nx*dg_els(ie)%nx + zt_coef*(1.d0 - dg_els(ie)%nx*dg_els(ie)%nx);
1027 compz(1,2) = zn_coef*dg_els(ie)%nx*dg_els(ie)%ny - zt_coef*dg_els(ie)%nx*dg_els(ie)%ny;
1028 compz(1,3) = zn_coef*dg_els(ie)%nx*dg_els(ie)%nz - zt_coef*dg_els(ie)%nx*dg_els(ie)%nz;
1029 compz(2,2) = zn_coef*dg_els(ie)%ny*dg_els(ie)%ny + zt_coef*(1.d0 - dg_els(ie)%ny*dg_els(ie)%ny);
1030 compz(2,3) = zn_coef*dg_els(ie)%ny*dg_els(ie)%nz - zt_coef*dg_els(ie)%ny*dg_els(ie)%nz;
1031 compz(3,3) = zn_coef*dg_els(ie)%nz*dg_els(ie)%nz + zt_coef*(1.d0 - dg_els(ie)%nz*dg_els(ie)%nz);
1032 compz(2,1) = compz(1,2); compz(3,1) = compz(1,3); compz(3,2) = compz(2,3);
1033
1034! write(*,*)
1035! write(*,*) compZ(1,1),compZ(1,2),compZ(1,3),compZ(2,1),compZ(2,2),compZ(2,3),compZ(3,1),compZ(3,2),compZ(3,3)
1036! read(*,*)
1037
1038 detz = compz(1,1)*compz(2,2)*compz(3,3) + compz(1,2)*compz(2,3)*compz(3,1) + compz(1,3)*compz(2,1)*compz(3,2) &
1039 - compz(1,3)*compz(2,2)*compz(3,1) - compz(1,1)*compz(2,3)*compz(3,2) - compz(1,2)*compz(2,1)*compz(3,3);
1040
1041! write(*,*) detZ;
1042! read(*,*)
1043
1044 z11 = compz(2,2)*compz(3,3) - compz(2,3)*compz(3,2); z11 = z11/detz;
1045 z12 = -compz(2,1)*compz(3,3) + compz(2,3)*compz(3,1); z12 = z12/detz;
1046 z13 = compz(2,1)*compz(2,3) - compz(2,2)*compz(3,1); z13 = z13/detz;
1047 z21 = -compz(1,2)*compz(3,3) + compz(1,3)*compz(3,2); z21 = z21/detz;
1048 z22 = compz(1,1)*compz(3,3) - compz(1,3)*compz(3,1); z22 = z22/detz;
1049 z23 = -compz(1,1)*compz(3,2) + compz(2,1)*compz(3,1); z23 = z23/detz;
1050 z31 = compz(1,2)*compz(2,3) - compz(1,3)*compz(2,2); z31 = z31/detz;
1051 z32 = -compz(1,1)*compz(2,3) + compz(1,3)*compz(2,1); z32 = z32/detz;
1052 z33 = compz(1,1)*compz(2,2) - compz(1,2)*compz(2,1); z33 = z33/detz;
1053
1054 invz(1,1) = z11; invz(1,2) = z21; invz(1,3) = z31;
1055 invz(2,1) = z12; invz(2,2) = z22; invz(2,3) = z32;
1056 invz(3,1) = z13; invz(3,2) = z23; invz(3,3) = z33;
1057
1058! Ide = matmul(invZ,compZ);
1059! write(*,*) Ide
1060! read(*,*)
1061
1062 endif
1063
1064
1065
1066 allocate(el_new(ie)%matP(3*(nn**3),3*(nn**3))); el_new(ie)%matP = 0.d0
1067 allocate(el_new(ie)%matM(dg_els(ie)%nofne))
1068 !--------------------------------------------------------------------------------------------
1069 if (testmode .eq. 1) then
1070 allocate(el_new(ie)%matP_only_uv(3*(nn**3),3*(nn**3))); el_new(ie)%matP_only_uv = 0.d0
1071 allocate(el_new(ie)%matM_only_uv(dg_els(ie)%nofne))
1072 endif
1073 !--------------------------------------------------------------------------------------------
1074
1075 do ic = 1, dg_els(ie)%nofne
1076
1077 imne = dg_els(ie)%conf(ic,0)
1078 iene = dg_els(ie)%conf(ic,1)
1079 ifacene = dg_els(ie)%conf(ic,2)
1080
1081 mm = 2
1082
1083 do i = 1, nm
1084 if(tag_mat(i) .eq. imne ) then
1085 mm = sd(i) +1
1086 lambda = 2.d0*prop_mat(im,2)*prop_mat(i,2)/(prop_mat(im,2) + prop_mat(i,2))
1087 mu = 2.d0*prop_mat(im,3)*prop_mat(i,3)/(prop_mat(im,3) + prop_mat(i,3))
1088 endif
1089 enddo
1090
1091 el_new(ie)%el_conf(ic,0) = imne
1092 el_new(ie)%el_conf(ic,1) = iene
1093 el_new(ie)%el_conf(ic,2) = ifacene
1094
1095 call get_face_dg(faces, nel_dg_glo, iene, ifacene, ic2)
1096
1097 allocate(jp(3*(nn**3),3*(nn**3)),el_new(ie)%matM(ic)%MJUMP(3*(nn**3),3*(mm**3)),jm(3*(nn**3),3*(mm**3)))
1098 !---------------------------------------------------------------------------------------------------
1099 if(testmode .eq. 1) then
1100 allocate(jp_only_uv(3*(nn**3),3*(nn**3)),&
1101 el_new(ie)%matM(ic)%MJUMP_only_uv(3*(nn**3),3*(mm**3)),jm_only_uv(3*(nn**3),3*(mm**3)))
1102 endif
1103 !---------------------------------------------------------------------------------------------------
1104
1105 pen_h = min(area_nodes(1,ic1), area_nodes(1,ic2))
1106 pen_p = max(nn-1,mm-1)
1107
1108 pen = penalty_c * (lambda + 2.d0*mu) * pen_p**2.d0 / pen_h
1109
1110 cp_a = 0.5d0*(lambda+2.d0*mu)*dg_els(ie)%nx
1111 cp_b = 0.5d0*lambda*dg_els(ie)%nx
1112 cp_c = 0.5d0*mu*dg_els(ie)%ny
1113 cp_d = 0.5d0*mu*dg_els(ie)%nz
1114 cp_e = 0.5d0*mu*dg_els(ie)%nx
1115 cp_f = 0.5d0*(lambda+2.d0*mu)*dg_els(ie)%ny
1116 cp_g = 0.5d0*lambda*dg_els(ie)%ny
1117 cp_n = 0.5d0*(lambda+2.d0*mu)*dg_els(ie)%nz
1118 cp_p = 0.5d0*lambda*dg_els(ie)%nz
1119
1120 call make_bilinear_map(area_nodes(2:25,ic2), &
1121 c_alfa11, c_alfa12, c_alfa13, &
1122 c_alfa21, c_alfa22, c_alfa23, &
1123 c_alfa31, c_alfa32, c_alfa33, &
1124 c_beta11, c_beta12, c_beta13, &
1125 c_beta21, c_beta22, c_beta23, &
1126 c_beta31, c_beta32, c_beta33, &
1127 c_gamma1, c_gamma2, c_gamma3, &
1128 c_delta1, c_delta2, c_delta3)
1129
1130
1131
1132 call make_loc_matrix_dg(dg_els(ie)%x_pl, dg_els(ie)%y_pl, dg_els(ie)%z_pl, &
1133 dg_els(ie)%wx_pl, dg_els(ie)%wy_pl, dg_els(ie)%wz_pl,&
1134 dg_els(ie)%x_mn, dg_els(ie)%y_mn, dg_els(ie)%z_mn, &
1135 dg_els(ie)%quad_rule, nn, mm, &
1136 dg_els(ie)%omega_minus(1:dg_els(ie)%quad_rule,0:3), &
1137 alfa11(ne1),alfa12(ne1),alfa13(ne1),&
1138 alfa21(ne1),alfa22(ne1),alfa23(ne1),&
1139 alfa31(ne1),alfa32(ne1),alfa33(ne1),&
1140 beta11(ne1),beta12(ne1),beta13(ne1),&
1141 beta21(ne1),beta22(ne1),beta23(ne1),&
1142 beta31(ne1),beta32(ne1),beta33(ne1),&
1143 gamma1(ne1),gamma2(ne1),gamma3(ne1),&
1144 delta1(ne1),delta2(ne1),delta3(ne1),&
1145 iene, c_alfa11, c_alfa12, c_alfa13,&
1146 c_alfa21,c_alfa22,c_alfa23,&
1147 c_alfa31,c_alfa32,c_alfa33,&
1148 c_beta11,c_beta12,c_beta13,&
1149 c_beta21,c_beta22,c_beta23,&
1150 c_beta31,c_beta32,c_beta33,&
1151 c_gamma1,c_gamma2,c_gamma3,&
1152 c_delta1,c_delta2,c_delta3,&
1153 cp_a,cp_b,cp_c,cp_d,cp_e,cp_f,cp_g,cp_n,cp_p, &
1154 pen, dg_cnst, jp, jm, mpi_id, pen_h, testmode,&
1155 jp_only_uv, jm_only_uv,fracture_yn,invz)
1156
1157 el_new(ie)%matP = el_new(ie)%matP + jp
1158 el_new(ie)%matM(ic)%MJUMP = jm
1159
1160 !--------------------------------------------------------------------
1161 if(testmode .eq. 1) then
1162 el_new(ie)%matP_only_uv = el_new(ie)%matP_only_uv + jp_only_uv
1163 el_new(ie)%matM(ic)%MJUMP_only_uv = jm_only_uv
1164 endif
1165 !--------------------------------------------------------------------
1166
1167
1168 do i = 1, 3*nn**3
1169 do j = 1, 3*mm**3
1170 if(el_new(ie)%matM(ic)%MJUMP(i,j) .ne. 0.d0) nnz_m = nnz_m + 1
1171 !------------------------------------------------------------------------------
1172 if(testmode .eq. 1) then
1173 if (el_new(ie)%matM(ic)%MJUMP_only_uv(i,j) .ne. 0.d0) &
1174 nnz_m_only_uv = nnz_m_only_uv + 1
1175 endif
1176 !------------------------------------------------------------------------------
1177 enddo
1178 enddo
1179
1180 el_new(ie)%nnz_col = el_new(ie)%nnz_col + 3*mm**3
1181
1182 !--------------------------------------------------------------------
1183 if(testmode .eq. 1) then
1184 el_new(ie)%nnz_col_only_uv = el_new(ie)%nnz_col_only_uv + 3*mm**3
1185 deallocate(jm_only_uv,jp_only_uv)
1186 endif
1187 !--------------------------------------------------------------------
1188
1189 deallocate(jm,jp)
1190 enddo
1191
1192
1193 do i = 1, 3*nn**3
1194 do j = 1, 3*nn**3
1195 if(el_new(ie)%matP(i,j) .ne. 0.d0) nnz_p = nnz_p + 1
1196 enddo
1197 enddo
1198
1199 el_new(ie)%nnz_minus = nnz_m; el_new(ie)%nnz_plus = nnz_p
1200
1201 !--------------------------------------------------------------------------------------------
1202 if(testmode .eq. 1) then
1203 do i = 1, 3*nn**3
1204 do j = 1, 3*nn**3
1205 if(el_new(ie)%matP_only_uv(i,j) .ne. 0.d0) nnz_p_only_uv = nnz_p_only_uv + 1
1206 enddo
1207 enddo
1208 el_new(ie)%nnz_minus_only_uv = nnz_m_only_uv; el_new(ie)%nnz_plus_only_uv = nnz_p_only_uv
1209 endif
1210 !---------------------------------------------------------------------------------------------
1211
1212!
1213! STORING MATRIX IN A SPARSE FORMAT
1214!
1215
1216 nn = el_new(ie)%deg
1217 allocate(el_new(ie)%IPlus(0:3*nn**3), el_new(ie)%IMin(0:3*nn**3))
1218 allocate(el_new(ie)%JPlus(el_new(ie)%nnz_plus),el_new(ie)%matPlus(el_new(ie)%nnz_plus))
1219 allocate(el_new(ie)%JMin(el_new(ie)%nnz_minus),el_new(ie)%matMin(el_new(ie)%nnz_minus))
1220 allocate(j4s(el_new(ie)%nnz_minus),m4s(el_new(ie)%nnz_minus))
1221 allocate(i4s(el_new(ie)%nnz_plus))
1222
1223 k = 1
1224 do i = 1, 3*nn**3
1225 do j = 1, 3*nn**3
1226 if( el_new(ie)%matP(i,j) .ne. 0.d0 ) then
1227 i4s(k) = i
1228 el_new(ie)%JPlus(k) = j
1229 el_new(ie)%matPlus(k) = el_new(ie)%matP(i,j)
1230 k = k + 1
1231 endif
1232 enddo
1233 enddo
1234
1235 deallocate(el_new(ie)%matP)
1236
1237 el_new(ie)%IPlus= 0
1238 do i = 1, el_new(ie)%nnz_plus
1239 el_new(ie)%IPlus(i4s(i)) = el_new(ie)%IPlus(i4s(i)) + 1
1240 enddo
1241 do i = 1, 3*nn**3
1242 el_new(ie)%IPlus(i) = el_new(ie)%IPlus(i) + el_new(ie)%IPlus(i-1)
1243 enddo
1244
1245 deallocate(i4s)
1246 allocate(i4s(el_new(ie)%nnz_minus))
1247
1248 k = 1
1249 ishift = 0
1250 jshift = 0
1251
1252 do ic = 1, el_new(ie)%num_of_ne
1253
1254 mm = 2
1255 do i = 1, nm
1256 if(tag_mat(i) .eq. el_new(ie)%el_conf(ic,0) ) then
1257 mm = sd(i) +1
1258 endif
1259 enddo
1260
1261
1262 do i = 1, 3*nn**3
1263 do j = 1, 3*mm**3
1264 if( el_new(ie)%matM(ic)%MJUMP(i,j) .ne. 0.d0 ) then
1265 i4s(k) = i
1266 j4s(k) = j + jshift
1267 m4s(k) = el_new(ie)%matM(ic)%MJUMP(i,j)
1268
1269 k = k + 1
1270 endif
1271
1272 enddo
1273 enddo
1274
1275
1276 jshift = jshift + 3*mm**3
1277
1278 deallocate(el_new(ie)%matM(ic)%MJUMP)
1279
1280
1281 enddo
1282
1283 k = 1
1284 do i = 1, 3*nn**3
1285 do j = 1, el_new(ie)%nnz_minus
1286 if(i4s(j) .eq. i) then
1287 el_new(ie)%JMin(k) = j4s(j)
1288 el_new(ie)%matMin(k) = m4s(j)
1289 k = k + 1
1290 endif
1291 enddo
1292 enddo
1293
1294
1295 deallocate(j4s, m4s)
1296
1297 el_new(ie)%IMin = 0
1298 do i = 1, el_new(ie)%nnz_minus
1299 el_new(ie)%IMin(i4s(i)) = el_new(ie)%IMin(i4s(i)) + 1
1300 enddo
1301
1302 do i = 1, 3*nn**3
1303 el_new(ie)%IMin(i) = el_new(ie)%IMin(i) + el_new(ie)%IMin(i-1)
1304 enddo
1305
1306 deallocate(i4s)
1307
1308 !-----------------------------------------------------------------------------------------
1309 if(testmode .eq. 1) then
1310
1311 nn = el_new(ie)%deg
1312 allocate(el_new(ie)%IPlus_only_uv(0:3*nn**3),el_new(ie)%IMin_only_uv(0:3*nn**3))
1313 allocate(el_new(ie)%JPlus_only_uv(el_new(ie)%nnz_plus_only_uv))
1314 allocate(el_new(ie)%matPlus_only_uv(el_new(ie)%nnz_plus_only_uv))
1315 allocate(el_new(ie)%JMin_only_uv(el_new(ie)%nnz_minus_only_uv))
1316 allocate(el_new(ie)%matMin_only_uv(el_new(ie)%nnz_minus_only_uv))
1317 allocate(j4s(el_new(ie)%nnz_minus_only_uv),m4s(el_new(ie)%nnz_minus_only_uv))
1318 allocate(i4s(el_new(ie)%nnz_plus_only_uv))
1319
1320 k = 1
1321 do i = 1, 3*nn**3
1322 do j = 1, 3*nn**3
1323 if( el_new(ie)%matP_only_uv(i,j) .ne. 0.d0 ) then
1324 i4s(k) = i
1325 el_new(ie)%JPlus_only_uv(k) = j
1326 el_new(ie)%matPlus_only_uv(k) = el_new(ie)%matP_only_uv(i,j)
1327 k = k + 1
1328 endif
1329 enddo
1330 enddo
1331
1332 deallocate(el_new(ie)%matP_only_uv)
1333
1334 el_new(ie)%IPlus_only_uv= 0
1335 do i = 1, el_new(ie)%nnz_plus_only_uv
1336 el_new(ie)%IPlus_only_uv(i4s(i)) = el_new(ie)%IPlus_only_uv(i4s(i)) + 1
1337 enddo
1338 do i = 1, 3*nn**3
1339 el_new(ie)%IPlus_only_uv(i) = el_new(ie)%IPlus_only_uv(i) + el_new(ie)%IPlus_only_uv(i-1)
1340 enddo
1341
1342 deallocate(i4s)
1343 allocate(i4s(el_new(ie)%nnz_minus_only_uv))
1344
1345 k = 1
1346 ishift = 0
1347 jshift = 0
1348
1349 do ic = 1, el_new(ie)%num_of_ne
1350
1351 mm = 2
1352 do i = 1, nm
1353 if(tag_mat(i) .eq. el_new(ie)%el_conf(ic,0) ) then
1354 mm = sd(i) +1
1355 endif
1356 enddo
1357
1358
1359 do i = 1, 3*nn**3
1360 do j = 1, 3*mm**3
1361 if( el_new(ie)%matM(ic)%MJUMP_only_uv(i,j) .ne. 0.d0 ) then
1362 i4s(k) = i
1363 j4s(k) = j + jshift
1364 m4s(k) = el_new(ie)%matM(ic)%MJUMP_only_uv(i,j)
1365
1366 k = k + 1
1367 endif
1368 enddo
1369 enddo
1370
1371 jshift = jshift + 3*mm**3
1372
1373 deallocate(el_new(ie)%matM(ic)%MJUMP_only_uv)
1374 enddo
1375
1376 k = 1
1377 do i = 1, 3*nn**3
1378 do j = 1, el_new(ie)%nnz_minus_only_uv
1379 if(i4s(j) .eq. i) then
1380 el_new(ie)%JMin_only_uv(k) = j4s(j)
1381 el_new(ie)%matMin_only_uv(k) = m4s(j)
1382 k = k + 1
1383 endif
1384 enddo
1385 enddo
1386
1387 deallocate(j4s, m4s)
1388
1389 el_new(ie)%IMin_only_uv = 0
1390 do i = 1, el_new(ie)%nnz_minus_only_uv
1391 el_new(ie)%IMin_only_uv(i4s(i)) = el_new(ie)%IMin_only_uv(i4s(i)) + 1
1392 enddo
1393
1394 do i = 1, 3*nn**3
1395 el_new(ie)%IMin_only_uv(i) = el_new(ie)%IMin_only_uv(i) + el_new(ie)%IMin_only_uv(i-1)
1396 enddo
1397
1398 deallocate(i4s)
1399
1400
1401 endif
1402 !-----------------------------------------------------------------------------------------
1403
1404 enddo
1405
1406!!$OMP END DO
1407!!$OMP END PARALLEL
1408
1409
1410 deallocate(dg_els)
1411 deallocate(con_dg_frc_int,con_dg_frc_real)
1412
1413
1414 return
1415
1416 end subroutine make_dg_interface
1417
1418
get_face_dg
subroutine get_face_dg(faces, nel_dg_glo, ik, ind, face_found)
Find the face (quad) of a DG element (hex)
Definition GET_FACE_DG.f90:34
get_indloc_from_indglo
subroutine get_indloc_from_indglo(local_el, nel_loc, ie, ic)
Returns local id from global id.
Definition GET_INDLOC_FROM_INDGLO.f90:30
get_neighbour_elem
subroutine get_neighbour_elem(omega_m, npt, n, id_mpi, n_qp)
Computes neighbouring element.
Definition GET_NEIGHBOUR_ELEM.f90:33
make_bilinear_map
subroutine make_bilinear_map(nodes, c_alfa11, c_alfa12, c_alfa13, c_alfa21, c_alfa22, c_alfa23, c_alfa31, c_alfa32, c_alfa33, c_beta11, c_beta12, c_beta13, c_beta21, c_beta22, c_beta23, c_beta31, c_beta32, c_beta33, c_gamma1, c_gamma2, c_gamma3, c_delta1, c_delta2, c_delta3)
Takes values from nodes array.
Definition MAKE_BILINEAR_MAP.f90:57
make_dg_interface
subroutine make_dg_interface(nm, sd, tag_mat, prop_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_cnst, penalty_c, faces, area_nodes, el_new, filename, testmode)
Creates data structure for dg elements.
Definition MAKE_DG_INTERFACE.f90:93
make_loc_matrix_dg
subroutine make_loc_matrix_dg(x_pl, y_pl, z_pl, wx_pl, wy_pl, wz_pl, x_mn, y_mn, z_mn, nq, nn, mm, o_minus, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, ine, alfa11_mn, alfa12_mn, alfa13_mn, alfa21_mn, alfa22_mn, alfa23_mn, alfa31_mn, alfa32_mn, alfa33_mn, beta11_mn, beta12_mn, beta13_mn, beta21_mn, beta22_mn, beta23_mn, beta31_mn, beta32_mn, beta33_mn, gamma1_mn, gamma2_mn, gamma3_mn, delta1_mn, delta2_mn, delta3_mn, cp_a, cp_b, cp_c, cp_d, cp_e, cp_f, cp_g, cp_n, cp_p, pen, dg_cnst, jp, jm, id_mpi, det_scal, testmode, jp_uv, jm_uv, fr_yn, invz)
Computes DG matrices for jumps.
Definition MAKE_LOC_MATRIX_DG.f90:145
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
newton_rapson
subroutine newton_rapson(xnod, ynod, znod, alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33, beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33, gamma1, gamma2, gamma3, delta1, delta2, delta3, tt, csi_s, eta_s, zeta_s, nofi, id_mpi, ipiu, imeno, toll1, toll2, flag)
Performs the Newton Rapson algorithm.
Definition NEWTON_RAPSON.f90:73
dgjump
Contains structure for jump matrices.
Definition MODULES.f90:155
max_var
Set maximal bounds.
Definition MODULES.f90:54
max_var::max_quad_points
integer, parameter max_quad_points
max number of quadrature nodes on a DG surface
Definition MODULES.f90:59
max_var::nofinr
integer, parameter nofinr
max number of newton rapson iterations
Definition MODULES.f90:58
str_mesh_after
Contains mesh structure for DG elements after pre-processing.
Definition MODULES.f90:115
dgjump::el4loop
element structure for time loop (RAM saving)
Definition MODULES.f90:167
str_mesh_after::element_after
Definition MODULES.f90:119