SPEED
MAKE_LOC_MATRIX_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
105
106
107 subroutine make_loc_matrix_dg(x_pl, y_pl, z_pl, &
108 wx_pl,wy_pl, wz_pl,&
109 x_mn, y_mn, z_mn, &
110 nq, nn, mm, &
111 o_minus, &
112 alfa11,alfa12,alfa13,&
113 alfa21,alfa22,alfa23,&
114 alfa31,alfa32,alfa33,&
115 beta11,beta12,beta13,&
116 beta21,beta22,beta23,&
117 beta31,beta32,beta33,&
118 gamma1,gamma2,gamma3,&
119 delta1,delta2,delta3,&
120 ine,&
121 alfa11_mn,alfa12_mn,alfa13_mn,&
122 alfa21_mn,alfa22_mn,alfa23_mn,&
123 alfa31_mn,alfa32_mn,alfa33_mn,&
124 beta11_mn,beta12_mn,beta13_mn,&
125 beta21_mn,beta22_mn,beta23_mn,&
126 beta31_mn,beta32_mn,beta33_mn,&
127 gamma1_mn,gamma2_mn,gamma3_mn,&
128 delta1_mn,delta2_mn,delta3_mn,&
129 cp_a,cp_b,cp_c,cp_d,cp_e,cp_f,cp_g,cp_n,cp_p, &
130 pen, dg_cnst,JP,JM,id_mpi,det_scal, testmode,&
131 JP_uv, JM_uv, fr_yn, invZ)
132
133
134 use max_var
135 use str_mesh
136
137 implicit none
138
139
140 integer*4 :: i, id_mpi, testmode, fr_yn
141 integer*4 :: ip1 , i1, j1, h1, m1, n1, p1, L1, K1
142 integer*4 :: ip2 , i2, j2, h2, m2, n2, p2, L2, K2
143 integer*4 :: ip3 , i3, j3, h3, m3, n3, p3, L3, K3
144 !integer*4 :: nthreads, OMP_GET_NUM_THREADS
145
146 integer*4, intent(in) :: nn, nq, mm, ine
147
148 integer*4, dimension(1:nq,0:3), intent(in) :: o_minus
149
150 real*8 :: phi, dphi
151 real*8 :: dxdx1,dxdy1,dxdz1,dydx1,dydy1,dydz1,dzdx1,dzdy1,dzdz1,det_j1, det_s1
152 real*8 :: dxdx2,dxdy2,dxdz2,dydx2,dydy2,dydz2,dzdx2,dzdy2,dzdz2,det_j2, det_s2
153 real*8 :: dxdx3,dxdy3,dxdz3,dydx3,dydy3,dydz3,dzdx3,dzdy3,dzdz3,det_j3, det_s3
154 real*8 :: dxdu1, dydu1, dzdu1, dxdv1, dydv1, dzdv1
155 real*8 :: dxdu2, dydu2, dzdu2, dxdv2, dydv2, dzdv2
156 real*8 :: dxdu3, dydu3, dzdu3, dxdv3, dydv3, dzdv3
157 real*8 :: dxdx_m,dxdy_m,dxdz_m,dydx_m,dydy_m,dydz_m,dzdx_m,dzdy_m,dzdz_m,det_j_m
158 real*8 :: dpdc1, dpde1, dpdz1, dpdx1, dpdy1, dpdzt1, plx1, pkx1
159 real*8 :: dpdc2, dpde2, dpdz2, dpdx2, dpdy2, dpdzt2, plx2, pkx2
160 real*8 :: dpdc3, dpde3, dpdz3, dpdx3, dpdy3, dpdzt3, plx3, pkx3
161 real*8 :: tempo1, tempo2
162
163 real*8, intent(in) :: alfa11, alfa12, alfa13, alfa21, alfa22, alfa23, alfa31, alfa32, alfa33
164 real*8, intent(in) :: beta11, beta12, beta13, beta21, beta22, beta23, beta31, beta32, beta33
165 real*8, intent(in) :: gamma1, gamma2, gamma3, delta1, delta2, delta3
166 real*8, intent(in) :: alfa11_mn, alfa12_mn, alfa13_mn
167 real*8, intent(in) :: alfa21_mn, alfa22_mn, alfa23_mn
168 real*8, intent(in) :: alfa31_mn, alfa32_mn, alfa33_mn
169 real*8, intent(in) :: beta11_mn, beta12_mn, beta13_mn
170 real*8, intent(in) :: beta21_mn, beta22_mn, beta23_mn
171 real*8, intent(in) :: beta31_mn, beta32_mn, beta33_mn
172 real*8, intent(in) :: gamma1_mn, gamma2_mn, gamma3_mn
173 real*8, intent(in) :: delta1_mn, delta2_mn, delta3_mn
174 real*8, intent(in) :: cp_a, cp_b, cp_c, cp_d, cp_e, cp_f, cp_g, cp_n, cp_p
175 real*8, intent(in) :: pen, dg_cnst, det_scal
176 real*8, dimension(nq), intent(in) :: x_pl, y_pl, z_pl, wx_pl, wy_pl, wz_pl, x_mn, y_mn, z_mn
177 real*8, dimension(nn) :: ctp
178 real*8, dimension(mm) :: ctm
179
180 real*8, dimension(:,:), allocatable :: dd,ee,ff,gg,hh,ii,qq
181 real*8, dimension(:,:), allocatable :: aa,bb,cc,pp
182 real*8, dimension(3*nn**3,3*nn**3), intent(out) :: jp, jp_uv
183 real*8, dimension(3*nn**3,3*mm**3), intent(out) :: jm, jm_uv
184 real*8, dimension(3,3) :: invz
185
186
187 call make_lgl_nodes(nn,ctp)
188 call make_lgl_nodes(mm,ctm)
189
190!********************************************************************************************!
191! MATRICES (+,+) !
192!********************************************************************************************!
193
194
195
196
197 allocate(aa(nn**3,nn**3),bb(nn**3,nn**3),cc(nn**3,nn**3),pp(nn**3,nn**3))
198
199 aa = 0.d0
200 bb = 0.d0
201 cc = 0.d0
202 pp = 0.d0
203
204
205 allocate(dd(nn**3,mm**3),ee(nn**3,mm**3),ff(nn**3,mm**3),qq(nn**3,mm**3))
206
207 dd = 0.d0
208 ee = 0.d0
209 ff = 0.d0
210 qq = 0.d0
211
212
213
214 allocate(gg(nn**3,mm**3),hh(nn**3,mm**3),ii(nn**3,mm**3))
215 gg = 0.d0
216 hh = 0.d0
217 ii = 0.d0
218
219
220
221!!$OMP PARALLEL
222!!$OMP SECTIONS
223
224!!$OMP SECTION
225
226
227 do p1 = 1, nn
228 do n1 = 1, nn
229 do m1 = 1, nn
230 l1 = m1 + (n1-1)*nn + (p1-1)*nn*nn
231
232
233 do h1 = 1, nn
234 do j1 = 1, nn
235 do i1 = 1, nn
236 k1 = i1 + (j1-1)*nn + (h1-1)*nn*nn
237
238 do ip1 = 1, nq
239
240 if( o_minus(ip1,1) .eq. ine) then
241
242 dxdx1 = alfa11 + beta12*z_pl(ip1) + beta13*y_pl(ip1) + gamma1*y_pl(ip1)*z_pl(ip1)
243 dydx1 = alfa21 + beta22*z_pl(ip1) + beta23*y_pl(ip1) + gamma2*y_pl(ip1)*z_pl(ip1)
244 dzdx1 = alfa31 + beta32*z_pl(ip1) + beta33*y_pl(ip1) + gamma3*y_pl(ip1)*z_pl(ip1)
245
246 dxdy1 = alfa12 + beta11*z_pl(ip1) + beta13*x_pl(ip1) + gamma1*z_pl(ip1)*x_pl(ip1)
247 dydy1 = alfa22 + beta21*z_pl(ip1) + beta23*x_pl(ip1) + gamma2*z_pl(ip1)*x_pl(ip1)
248 dzdy1 = alfa32 + beta31*z_pl(ip1) + beta33*x_pl(ip1) + gamma3*z_pl(ip1)*x_pl(ip1)
249
250 dxdz1 = alfa13 + beta11*y_pl(ip1) + beta12*x_pl(ip1) + gamma1*x_pl(ip1)*y_pl(ip1)
251 dydz1 = alfa23 + beta21*y_pl(ip1) + beta22*x_pl(ip1) + gamma2*x_pl(ip1)*y_pl(ip1)
252 dzdz1 = alfa33 + beta31*y_pl(ip1) + beta32*x_pl(ip1) + gamma3*x_pl(ip1)*y_pl(ip1)
253
254 det_j1 = dxdz1 * (dydx1*dzdy1 - dzdx1*dydy1) &
255 - dydz1 * (dxdx1*dzdy1 - dzdx1*dxdy1) &
256 + dzdz1 * (dxdx1*dydy1 - dydx1*dxdy1)
257
258
259 dpdc1 = dphi(nn-1, ctp(i1), x_pl(ip1)) * phi(nn-1, ctp(j1), y_pl(ip1)) * phi(nn-1, ctp(h1), z_pl(ip1))
260 dpde1 = phi(nn-1, ctp(i1), x_pl(ip1)) * dphi(nn-1, ctp(j1), y_pl(ip1)) * phi(nn-1, ctp(h1), z_pl(ip1))
261 dpdz1 = phi(nn-1, ctp(i1), x_pl(ip1)) * phi(nn-1, ctp(j1), y_pl(ip1)) * dphi(nn-1, ctp(h1), z_pl(ip1))
262
263 plx1 = phi(nn-1, ctp(m1), x_pl(ip1)) * phi(nn-1, ctp(n1), y_pl(ip1)) * phi(nn-1, ctp(p1), z_pl(ip1))
264 pkx1 = phi(nn-1, ctp(i1), x_pl(ip1)) * phi(nn-1, ctp(j1), y_pl(ip1)) * phi(nn-1, ctp(h1), z_pl(ip1))
265
266
267 dpdx1 = (dydx1*(dzdy1*dpdz1 - dzdz1*dpde1) - dydy1*(dzdx1*dpdz1 - dzdz1*dpdc1) + &
268 dydz1*(dzdx1*dpde1 - dzdy1*dpdc1))/det_j1
269
270 dpdy1 = (dzdx1*(dxdy1*dpdz1 - dxdz1*dpde1) - dzdy1*(dxdx1*dpdz1 - dxdz1*dpdc1) + &
271 dzdz1*(dxdx1*dpde1 - dxdy1*dpdc1))/det_j1
272
273
274 dpdzt1 = (dxdx1*(dydy1*dpdz1 - dydz1*dpde1) - dxdy1*(dydx1*dpdz1 - dydz1*dpdc1) + &
275 dxdz1*(dydx1*dpde1 - dydy1*dpdc1))/det_j1
276
277
278 det_s1 = det_scal/4.d0
279
280
281 aa(l1,k1) = aa(l1,k1) + det_s1 * wx_pl(ip1)*wy_pl(ip1)*wz_pl(ip1) * dpdx1 * plx1
282 bb(l1,k1) = bb(l1,k1) + det_s1 * wx_pl(ip1)*wy_pl(ip1)*wz_pl(ip1) * dpdy1 * plx1
283 cc(l1,k1) = cc(l1,k1) + det_s1 * wx_pl(ip1)*wy_pl(ip1)*wz_pl(ip1) * dpdzt1 * plx1
284 pp(l1,k1) = pp(l1,k1) + det_s1 * wx_pl(ip1)*wy_pl(ip1)*wz_pl(ip1) * pkx1 * plx1
285
286
287
288 endif
289
290
291 enddo ! end loop on ip1
292
293 enddo
294 enddo
295 enddo !end loop on K
296
297
298 enddo
299 enddo
300 enddo !end loop in L
301
302!********************************************************************************************!
303! MATRICES (-,+) !
304!********************************************************************************************!
305
306!!$OMP SECTION
307
308 do p2 = 1, nn
309 do n2 = 1, nn
310 do m2 = 1, nn
311 l2 = m2 + (n2-1)*nn + (p2-1)*nn*nn
312
313
314 do h2 = 1, mm
315 do j2 = 1, mm
316 do i2 = 1, mm
317 k2 = i2 + (j2-1)*mm + (h2-1)*mm*mm
318
319
320
321
322 do ip2 = 1, nq
323
324 if( o_minus(ip2,1) .eq. ine) then
325
326 dxdx2 = alfa11 + beta12*z_pl(ip2) + beta13*y_pl(ip2) + gamma1*y_pl(ip2)*z_pl(ip2)
327 dydx2 = alfa21 + beta22*z_pl(ip2) + beta23*y_pl(ip2) + gamma2*y_pl(ip2)*z_pl(ip2)
328 dzdx2 = alfa31 + beta32*z_pl(ip2) + beta33*y_pl(ip2) + gamma3*y_pl(ip2)*z_pl(ip2)
329
330 dxdy2 = alfa12 + beta11*z_pl(ip2) + beta13*x_pl(ip2) + gamma1*z_pl(ip2)*x_pl(ip2)
331 dydy2 = alfa22 + beta21*z_pl(ip2) + beta23*x_pl(ip2) + gamma2*z_pl(ip2)*x_pl(ip2)
332 dzdy2 = alfa32 + beta31*z_pl(ip2) + beta33*x_pl(ip2) + gamma3*z_pl(ip2)*x_pl(ip2)
333
334 dxdz2 = alfa13 + beta11*y_pl(ip2) + beta12*x_pl(ip2) + gamma1*x_pl(ip2)*y_pl(ip2)
335 dydz2 = alfa23 + beta21*y_pl(ip2) + beta22*x_pl(ip2) + gamma2*x_pl(ip2)*y_pl(ip2)
336 dzdz2 = alfa33 + beta31*y_pl(ip2) + beta32*x_pl(ip2) + gamma3*x_pl(ip2)*y_pl(ip2)
337
338
339 dxdx_m = alfa11_mn + beta12_mn*z_mn(ip2) + beta13_mn*y_mn(ip2) + gamma1_mn*y_mn(ip2)*z_mn(ip2)
340 dydx_m = alfa21_mn + beta22_mn*z_mn(ip2) + beta23_mn*y_mn(ip2) + gamma2_mn*y_mn(ip2)*z_mn(ip2)
341 dzdx_m = alfa31_mn + beta32_mn*z_mn(ip2) + beta33_mn*y_mn(ip2) + gamma3_mn*y_mn(ip2)*z_mn(ip2)
342
343 dxdy_m = alfa12_mn + beta11_mn*z_mn(ip2) + beta13_mn*x_mn(ip2) + gamma1_mn*z_mn(ip2)*x_mn(ip2)
344 dydy_m = alfa22_mn + beta21_mn*z_mn(ip2) + beta23_mn*x_mn(ip2) + gamma2_mn*z_mn(ip2)*x_mn(ip2)
345 dzdy_m = alfa32_mn + beta31_mn*z_mn(ip2) + beta33_mn*x_mn(ip2) + gamma3_mn*z_mn(ip2)*x_mn(ip2)
346
347 dxdz_m = alfa13_mn + beta11_mn*y_mn(ip2) + beta12_mn*x_mn(ip2) + gamma1_mn*x_mn(ip2)*y_mn(ip2)
348 dydz_m = alfa23_mn + beta21_mn*y_mn(ip2) + beta22_mn*x_mn(ip2) + gamma2_mn*x_mn(ip2)*y_mn(ip2)
349 dzdz_m = alfa33_mn + beta31_mn*y_mn(ip2) + beta32_mn*x_mn(ip2) + gamma3_mn*x_mn(ip2)*y_mn(ip2)
350
351 det_j_m = dxdz_m * (dydx_m*dzdy_m - dzdx_m*dydy_m) &
352 - dydz_m * (dxdx_m*dzdy_m - dzdx_m*dxdy_m) &
353 + dzdz_m * (dxdx_m*dydy_m - dydx_m*dxdy_m)
354
355
356
357 dpdc2 = dphi(mm-1, ctm(i2), x_mn(ip2)) * phi(mm-1, ctm(j2), y_mn(ip2)) * phi(mm-1, ctm(h2), z_mn(ip2))
358 dpde2 = phi(mm-1, ctm(i2), x_mn(ip2)) * dphi(mm-1, ctm(j2), y_mn(ip2)) * phi(mm-1, ctm(h2), z_mn(ip2))
359 dpdz2 = phi(mm-1, ctm(i2), x_mn(ip2)) * phi(mm-1, ctm(j2), y_mn(ip2)) * dphi(mm-1, ctm(h2), z_mn(ip2))
360
361
362 plx2 = phi(nn-1, ctp(m2), x_pl(ip2)) * phi(nn-1, ctp(n2), y_pl(ip2)) * phi(nn-1, ctp(p2), z_pl(ip2))
363 pkx2 = phi(mm-1, ctm(i2), x_mn(ip2)) * phi(mm-1, ctm(j2), y_mn(ip2)) * phi(mm-1, ctm(h2), z_mn(ip2))
364
365
366 dpdx2 = (dydx_m*(dzdy_m*dpdz2 - dzdz_m*dpde2) - dydy_m*(dzdx_m*dpdz2 - dzdz_m*dpdc2) + &
367 dydz_m*(dzdx_m*dpde2 - dzdy_m*dpdc2))/det_j_m
368
369 dpdy2 = (dzdx_m*(dxdy_m*dpdz2 - dxdz_m*dpde2) - dzdy_m*(dxdx_m*dpdz2 - dxdz_m*dpdc2) + &
370 dzdz_m*(dxdx_m*dpde2 - dxdy_m*dpdc2))/det_j_m
371
372 dpdzt2 = (dxdx_m*(dydy_m*dpdz2 - dydz_m*dpde2) - dxdy_m*(dydx_m*dpdz2 - dydz_m*dpdc2) + &
373 dxdz_m*(dydx_m*dpde2 - dydy_m*dpdc2))/det_j_m
374
375
376 det_s2 = det_scal/4.d0
377
378 dd(l2,k2) = dd(l2,k2) + det_s2 * wx_pl(ip2)*wy_pl(ip2)*wz_pl(ip2) * dpdx2 * plx2
379 ee(l2,k2) = ee(l2,k2) + det_s2 * wx_pl(ip2)*wy_pl(ip2)*wz_pl(ip2) * dpdy2 * plx2
380 ff(l2,k2) = ff(l2,k2) + det_s2 * wx_pl(ip2)*wy_pl(ip2)*wz_pl(ip2) * dpdzt2 * plx2
381 qq(l2,k2) = qq(l2,k2) + det_s2 * wx_pl(ip2)*wy_pl(ip2)*wz_pl(ip2) * pkx2 * plx2
382
383
384
385 endif
386
387
388 enddo ! end loop on ip2
389
390 enddo
391 enddo
392 enddo !end loop on K
393
394 enddo
395 enddo
396 enddo !end loop on L
397
398
399!********************************************************************************************!
400! MATRICES (+,-) !
401!********************************************************************************************!
402
403!!$OMP SECTION
404
405 do p3 = 1, nn
406 do n3 = 1, nn
407 do m3 = 1, nn
408 l3 = m3 + (n3-1)*nn + (p3-1)*nn*nn
409
410
411 do h3 = 1, mm
412 do j3 = 1, mm
413 do i3 = 1, mm
414 k3 = i3 + (j3-1)*mm + (h3-1)*mm*mm
415
416 do ip3 = 1, nq
417
418 if( o_minus(ip3,1) .eq. ine) then
419
420 dxdx3 = alfa11 + beta12*z_pl(ip3) + beta13*y_pl(ip3) + gamma1*y_pl(ip3)*z_pl(ip3)
421 dydx3 = alfa21 + beta22*z_pl(ip3) + beta23*y_pl(ip3) + gamma2*y_pl(ip3)*z_pl(ip3)
422 dzdx3 = alfa31 + beta32*z_pl(ip3) + beta33*y_pl(ip3) + gamma3*y_pl(ip3)*z_pl(ip3)
423
424 dxdy3 = alfa12 + beta11*z_pl(ip3) + beta13*x_pl(ip3) + gamma1*z_pl(ip3)*x_pl(ip3)
425 dydy3 = alfa22 + beta21*z_pl(ip3) + beta23*x_pl(ip3) + gamma2*z_pl(ip3)*x_pl(ip3)
426 dzdy3 = alfa32 + beta31*z_pl(ip3) + beta33*x_pl(ip3) + gamma3*z_pl(ip3)*x_pl(ip3)
427
428 dxdz3 = alfa13 + beta11*y_pl(ip3) + beta12*x_pl(ip3) + gamma1*x_pl(ip3)*y_pl(ip3)
429 dydz3 = alfa23 + beta21*y_pl(ip3) + beta22*x_pl(ip3) + gamma2*x_pl(ip3)*y_pl(ip3)
430 dzdz3 = alfa33 + beta31*y_pl(ip3) + beta32*x_pl(ip3) + gamma3*x_pl(ip3)*y_pl(ip3)
431
432 det_j3 = dxdz3 * (dydx3*dzdy3 - dzdx3*dydy3) &
433 - dydz3 * (dxdx3*dzdy3 - dzdx3*dxdy3) &
434 + dzdz3 * (dxdx3*dydy3 - dydx3*dxdy3)
435
436
437 dpdc3 = dphi(nn-1, ctp(m3), x_pl(ip3)) * phi(nn-1, ctp(n3), y_pl(ip3)) * phi(nn-1, ctp(p3), z_pl(ip3))
438 dpde3 = phi(nn-1, ctp(m3), x_pl(ip3)) * dphi(nn-1, ctp(n3), y_pl(ip3)) * phi(nn-1, ctp(p3), z_pl(ip3))
439 dpdz3 = phi(nn-1, ctp(m3), x_pl(ip3)) * phi(nn-1, ctp(n3), y_pl(ip3)) * dphi(nn-1, ctp(p3), z_pl(ip3))
440
441
442
443 pkx3 = phi(mm-1, ctm(i3), x_mn(ip3)) * phi(mm-1, ctm(j3), y_mn(ip3)) * phi(mm-1, ctm(h3), z_mn(ip3))
444
445
446 dpdx3 = (dydx3*(dzdy3*dpdz3 - dzdz3*dpde3) - dydy3*(dzdx3*dpdz3 - dzdz3*dpdc3) + &
447 dydz3*(dzdx3*dpde3 - dzdy3*dpdc3))/det_j3
448
449 dpdy3 = (dzdx3*(dxdy3*dpdz3 - dxdz3*dpde3) - dzdy3*(dxdx3*dpdz3 - dxdz3*dpdc3) + &
450 dzdz3*(dxdx3*dpde3 - dxdy3*dpdc3))/det_j3
451
452
453 dpdzt3 = (dxdx3*(dydy3*dpdz3 - dydz3*dpde3) - dxdy3*(dydx3*dpdz3 - dydz3*dpdc3) + &
454 dxdz3*(dydx3*dpde3 - dydy3*dpdc3))/det_j3
455
456
457
458 det_s3 = det_scal/4.d0
459
460 gg(l3,k3) = gg(l3,k3) + det_s3 * wx_pl(ip3)*wy_pl(ip3)*wz_pl(ip3) * dpdx3 * pkx3
461 hh(l3,k3) = hh(l3,k3) + det_s3 * wx_pl(ip3)*wy_pl(ip3)*wz_pl(ip3) * dpdy3 * pkx3
462 ii(l3,k3) = ii(l3,k3) + det_s3 * wx_pl(ip3)*wy_pl(ip3)*wz_pl(ip3) * dpdzt3 * pkx3
463
464
465
466
467 endif
468
469
470 enddo ! end loop on ip3
471
472 enddo
473 enddo
474 enddo !end loop on K
475
476 enddo
477 enddo
478 enddo !end loop on L
479
480
481!!$OMP END SECTIONS
482!!$OMP BARRIER
483!!$OMP END PARALLEL
484 if (fr_yn .eq. 0) then
485
486 jp = 0.d0
487 jp(1 : nn**3, 1 : nn**3) = pen*pp - cp_a*aa - cp_c*bb - cp_d*cc &
488 + dg_cnst*cp_a*transpose(aa) + dg_cnst*cp_c*transpose(bb) + dg_cnst*cp_d*transpose(cc)
489
490 jp(1 : nn**3, nn**3+1 : 2*nn**3) = - cp_c*aa - cp_b*bb &
491 + dg_cnst*cp_g*transpose(aa) + dg_cnst*cp_e*transpose(bb)
492
493 jp(1 : nn**3, 2*nn**3+1 : 3*nn**3) = - cp_d*aa - cp_b*cc &
494 + dg_cnst*cp_p*transpose(aa) + dg_cnst*cp_e*transpose(cc)
495
496! write(*,*) JP(1 : nn**3, 2*nn**3+1 : 3*nn**3)
497
498
499 jp(nn**3+1 : 2*nn**3, 1 : nn**3) = - cp_g*aa - cp_e*bb &
500 + dg_cnst*cp_c*transpose(aa) + dg_cnst*cp_b*transpose(bb)
501
502 jp(nn**3+1 : 2*nn**3, nn**3+1 : 2*nn**3) = pen*pp - cp_e*aa - cp_f*bb - cp_d*cc &
503 + dg_cnst*cp_e*transpose(aa) + dg_cnst*cp_f*transpose(bb) + dg_cnst*cp_d*transpose(cc)
504
505 jp(nn**3+1 : 2*nn**3, 2*nn**3+1 : 3*nn**3) = - cp_d*bb - cp_g*cc &
506 + dg_cnst*cp_p*transpose(bb) + dg_cnst*cp_c*transpose(cc)
507
508 jp(2*nn**3+1 : 3*nn**3, 1 : nn**3) = - cp_p*aa - cp_e*cc &
509 + dg_cnst*cp_d*transpose(aa) + dg_cnst*cp_b*transpose(cc)
510
511 jp(2*nn**3+1 : 3*nn**3, nn**3+1 : 2*nn**3) = - cp_p*bb - cp_c*cc &
512 + dg_cnst*cp_d*transpose(bb) + dg_cnst*cp_g*transpose(cc)
513
514 jp(2*nn**3+1 : 3*nn**3, 2*nn**3+1 : 3*nn**3) = pen*pp - cp_e*aa - cp_c*bb - cp_n*cc &
515 + dg_cnst*cp_e*transpose(aa) + dg_cnst*cp_c*transpose(bb) + dg_cnst*cp_n*transpose(cc)
516
517
518
519
520 jm = 0.d0
521 jm(1 : nn**3, 1 : mm**3) = - pen*qq - cp_a*dd - cp_c*ee - cp_d*ff &
522 - dg_cnst*cp_a*gg -dg_cnst*cp_c*hh - dg_cnst*cp_d*ii
523
524 jm(1 : nn**3, mm**3+1 : 2*mm**3) = - cp_c*dd - cp_b*ee - dg_cnst*cp_g*gg - dg_cnst*cp_e*hh
525
526 jm(1 : nn**3, 2*mm**3+1 : 3*mm**3) = - cp_d*dd - cp_b*ff - dg_cnst*cp_p*gg - dg_cnst*cp_e*ii
527
528 jm(nn**3+1 : 2*nn**3, 1 : mm**3) = - cp_g*dd - cp_e*ee - dg_cnst*cp_c*gg - dg_cnst*cp_b*hh
529
530 jm(nn**3+1 : 2*nn**3, mm**3+1 : 2*mm**3) = - pen*qq - cp_e*dd - cp_f*ee - cp_d*ff &
531 - dg_cnst*cp_e*gg - dg_cnst*cp_f*hh - dg_cnst*cp_d*ii
532
533 jm(nn**3+1 : 2*nn**3, 2*mm**3+1 : 3*mm**3) = - cp_d*ee - cp_g*ff - dg_cnst*cp_p*hh - dg_cnst*cp_c*ii
534
535 jm(2*nn**3+1 : 3*nn**3, 1 : mm**3) = - cp_p*dd - cp_e*ff - dg_cnst*cp_d*gg - dg_cnst*cp_b*ii
536
537 jm(2*nn**3+1 : 3*nn**3, mm**3+1 : 2*mm**3) = - cp_p*ee - cp_c*ff - dg_cnst*cp_d*hh - dg_cnst*cp_g*ii
538
539 jm(2*nn**3+1 : 3*nn**3, 2*mm**3+1 : 3*mm**3) = - pen*qq - cp_e*dd - cp_c*ee - cp_n*ff &
540 - dg_cnst*cp_e*gg - dg_cnst*cp_c*hh - dg_cnst*cp_n*ii
541
542 else
543
544 jp = 0.d0
545 jp(1 : nn**3, 1 : nn**3) = invz(1,1)*pp
546 jp(1 : nn**3, nn**3+1 : 2*nn**3) = invz(1,2)*pp
547 jp(1 : nn**3, 2*nn**3+1 : 3*nn**3) = invz(1,3)*pp
548 jp(nn**3+1 : 2*nn**3, 1 : nn**3) = invz(2,1)*pp
549 jp(nn**3+1 : 2*nn**3, nn**3+1 : 2*nn**3) = invz(2,2)*pp
550 jp(nn**3+1 : 2*nn**3, 2*nn**3+1 : 3*nn**3) = invz(2,3)*pp
551 jp(2*nn**3+1 : 3*nn**3, 1 : nn**3) = invz(3,1)*pp
552 jp(2*nn**3+1 : 3*nn**3, nn**3+1 : 2*nn**3) = invz(3,2)*pp
553 jp(2*nn**3+1 : 3*nn**3, 2*nn**3+1 : 3*nn**3) = invz(3,3)*pp
554
555
556
557 jm = 0.d0
558 jm(1 : nn**3, 1 : mm**3) = - invz(1,1)*qq
559 jm(1 : nn**3, mm**3+1 : 2*mm**3) = - invz(1,2)*qq
560 jm(1 : nn**3, 2*mm**3+1 : 3*mm**3) = - invz(1,3)*qq
561 jm(nn**3+1 : 2*nn**3, 1 : mm**3) = - invz(2,1)*qq
562 jm(nn**3+1 : 2*nn**3, mm**3+1 : 2*mm**3) = - invz(2,2)*qq
563 jm(nn**3+1 : 2*nn**3, 2*mm**3+1 : 3*mm**3) = - invz(2,3)*qq
564 jm(2*nn**3+1 : 3*nn**3, 1 : mm**3) = - invz(3,1)*qq
565 jm(2*nn**3+1 : 3*nn**3, mm**3+1 : 2*mm**3) = - invz(3,2)*qq
566 jm(2*nn**3+1 : 3*nn**3, 2*mm**3+1 : 3*mm**3) = - invz(3,3)*qq
567
568
569endif
570
571
572
573 if(testmode .eq. 1) then
574
575 jp_uv = 0.d0
576 jp_uv(1 : nn**3, 1 : nn**3) = pen*pp
577 jp_uv(nn**3+1 : 2*nn**3, nn**3+1 : 2*nn**3) = pen*pp
578 jp_uv(2*nn**3+1 : 3*nn**3, 2*nn**3+1 : 3*nn**3) = pen*pp
579
580 jm_uv = 0.d0
581 jm_uv(1 : nn**3, 1 : mm**3) = - pen*qq
582 jm_uv(nn**3+1 : 2*nn**3, mm**3+1 : 2*mm**3) = - pen*qq
583 jm_uv(2*nn**3+1 : 3*nn**3, 2*mm**3+1 : 3*mm**3) = - pen*qq
584
585 endif
586
587
588
589
590
591 deallocate(aa,bb,cc,pp)
592 deallocate(dd,ee,ff,qq)
593 deallocate(gg,hh,ii)
594
595
596 end subroutine make_loc_matrix_dg
make_lgl_nodes
subroutine make_lgl_nodes(np, ct)
Makes Gauss-Legendre-Lobatto nodes.
Definition MAKE_LGL_NODES.f90:27
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
phi
real *8 function phi(ng, csii, csij)
Evaluates Lengendre basis function on a specific point.
Definition PHI.f90:31
max_var
Set maximal bounds.
Definition MODULES.f90:54
str_mesh
Contains mesh structure for DG interface elements.
Definition MODULES.f90:84