Computes external loads.
96 nm,tag_mat,type_mat,sdeg_mat,tref_mat,prop_mat
97 ne_loc,local_el_num, &
98 alfa11,alfa12,alfa13,alfa21,alfa22,alfa23,&
99 alfa31,alfa32,alfa33,beta11,beta12,beta13,&
100 beta21,beta22,beta23,beta31,beta32,beta33,&
101 gamma1,gamma2,gamma3,delta1,delta2,delta3,&
102 cs_nnz_bc_loc,cs_bc_loc,&
103 nl_dirx,val_dirx,fun_dirx,tag_dirx,&
104 nl_diry,val_diry,fun_diry,tag_diry,&
105 nl_dirz,val_dirz,fun_dirz,tag_dirz,&
106 nl_neux,val_neux,fun_neux,tag_neux,&
107 nl_neuy,val_neuy,fun_neuy,tag_neuy,&
108 nl_neuz,val_neuz,fun_neuz,tag_neuz,&
109 nl_neun,val_neun,fun_neun,tag_neun,&
110 nl_poix,val_poix,fun_poix,&
111 nl_poiy,val_poiy,fun_poiy,&
112 nl_poiz,val_poiz,fun_poiz,&
113 nl_trax,val_trax,fun_trax,&
114 nl_tray,val_tray,fun_tray,&
115 nl_traz,val_traz,fun_traz,&
116 nl_plax,val_plax,fun_plax,tag_plax,&
117 nl_play,val_play,fun_play,tag_play,&
118 nl_plaz,val_plaz,fun_plaz,tag_plaz,&
119 srcmodflag,szsism, &
120 nl_sism,val_sism,fun_sism,tag_sism,&
121 nl_expl,val_expl,fun_expl,tag_expl,&
122 nl_forx,val_forx,fun_forx,&
123 nl_fory,val_fory,fun_fory,&
124 nl_forz,val_forz,fun_forz,&
125 nl_forc,val_forc,fun_forc,&
126 nl_pres,val_pres,fun_pres,&
127 nl_shea,val_shea,fun_shea,&
128 n_test,fun_test, &
129 nfunc,tag_func,func_type,func_indx,func_data,nfunc_data
130 fmat,&
131 con_hexa, nhexa,&
132 length_cns,&
133 sour_ns,max_num_ns,num_ns,&
134 facsmom,node_index_seq,&
135 tausmom,&
136 length_cne,&
137 sour_ne,max_num_ne,num_ne,&
138 facsexpl,&
139 mpi_comm, mpi_np, mpi_id,testmode, &
140 nmat_nhe, rho_nhe, lambda_nhe, mu_nhe)
141
143
144 implicit none
145
146 include 'SPEED.MPI'
147
148 character*12 :: name_prop
149 character*70 :: filename
150
151 integer*4 :: nnod_loc,cs_nnz_loc,nm,ne_loc,cs_nnz_bc_loc
152 integer*4 :: nl_dirX,nl_dirY,nl_dirZ,nl_neuX,nl_neuY,nl_neuZ
153 integer*4 :: nl_neuN
154integer*4 :: nl_poiX,nl_poiY,nl_poiZ,nl_forX,nl_forY,nl_forZ
155 integer*4 :: nl_plaX,nl_plaY,nl_plaZ,nl_traX,nl_traY,nl_traZ
156integer*4 :: nl_sism, srcmodflag, szsism
157 integer*4 :: nl_expl
158integer*4 :: nl_forc,nl_pres,nl_shea, n_test
159 integer*4 :: nfunc, nfunc_data
160 integer*4 :: im,ifun,ie,ip,il,nface_loc,nn,fn
161 integer*4 :: is,in,id
162 integer*4 :: i,j,k, val_st
163 integer*4 :: ipl, nb_tra_load
164integer*4 :: nhexa
165integer*4 :: sum_node_bottom,sum_node_bottom1,sum_node_bottom2
166integer*4 :: sum_node_bottom_first3
167integer*4 :: last_node_bottom,last_node_bottom1,last_node_bottom2
168integer*4 :: C,sit
169integer*4 :: isism,ipsism
170integer*4 :: length_cns, sum_length_cns
171integer*4 :: max_num_ns
172integer*4 :: tt, ie_glob, ie_surf
173 integer*4 :: ic, ic1, ic2, ic3, ic4
174 integer*4 :: mpi_comm, mpi_np, mpi_ierr, mpi_id
175 integer*4 :: iexpl,ipexpl
176 integer*4 :: length_cne, sum_length_cne
177 integer*4 :: max_num_ne
178integer*4 :: nnode_neuN,nelem_neuN
179 integer*4 :: ne1,ne2,ne3,ne4
180integer*4 :: index_vector,index, testmode, prova, nmat_nhe
181
182 integer*4, dimension(:), allocatable :: i4count
183 integer*4, dimension(:), allocatable :: i4normal
184 integer*4, dimension(:), allocatable :: ind_locX, ind_locY, ind_locZ
185 integer*4, dimension(:), allocatable :: ind_gloX, ind_gloY, ind_gloZ
186 integer*4, dimension(:), allocatable :: n_tra, n_tra_glo
187 integer*4, dimension(nnod_loc) :: node_index_seq
188integer*4, dimension(nl_expl) :: num_ne
189integer*4, dimension(6) :: num_sit
190integer*4, dimension(nl_sism) :: num_ns
191integer*4, dimension(nnod_loc) :: local_n_num
192 integer*4, dimension(ne_loc) :: local_el_num
193 integer*4, dimension(0:cs_nnz_loc) :: cs_loc
194 integer*4, dimension(nm) :: tag_mat,type_mat,sdeg_mat
195 integer*4, dimension(0:cs_nnz_bc_loc) :: cs_bc_loc
196 integer*4, dimension(nl_dirX) :: fun_dirX, tag_dirX
197 integer*4, dimension(nl_dirY) :: fun_dirY, tag_dirY
198 integer*4, dimension(nl_dirZ) :: fun_dirZ, tag_dirZ
199 integer*4, dimension(nl_neuX) :: fun_neuX, tag_neuX
200 integer*4, dimension(nl_neuY) :: fun_neuY, tag_neuY
201 integer*4, dimension(nl_neuZ) :: fun_neuZ, tag_neuZ
202 integer*4, dimension(nl_neuN) :: fun_neuN, tag_neuN
203 integer*4, dimension(nl_poiX) :: fun_poiX
204 integer*4, dimension(nl_poiY) :: fun_poiY
205 integer*4, dimension(nl_poiZ) :: fun_poiZ
206 integer*4, dimension(nl_traX) :: fun_traX
207 integer*4, dimension(nl_traY) :: fun_traY
208 integer*4, dimension(nl_traZ) :: fun_traZ
209 integer*4, dimension(nl_plaX) :: fun_plaX
210 integer*4, dimension(nl_plaY) :: fun_plaY
211 integer*4, dimension(nl_plaZ) :: fun_plaZ
212 integer*4, dimension(nl_plaX) :: tag_plaX
213integer*4, dimension(nl_plaY) :: tag_plaY
214 integer*4, dimension(nl_plaZ) :: tag_plaZ
215 integer*4, dimension(nl_sism) :: fun_sism
216 integer*4, dimension(nl_sism) :: tag_sism
217 integer*4, dimension(nl_expl) :: fun_expl
218 integer*4, dimension(nl_expl) :: tag_expl
219 integer*4, dimension(nl_forX) :: fun_forX
220 integer*4, dimension(nl_forY) :: fun_forY
221 integer*4, dimension(nl_forZ) :: fun_forZ
222 integer*4, dimension(nl_forc) :: fun_forc
223 integer*4, dimension(nl_pres) :: fun_pres
224 integer*4, dimension(n_test) :: fun_test
225 integer*4, dimension(nl_shea) :: fun_shea
226 integer*4, dimension(nfunc) :: tag_func
227 integer*4, dimension(nl_poiX) :: node_poiX
228 integer*4, dimension(nl_poiY) :: node_poiY
229 integer*4, dimension(nl_poiZ) :: node_poiZ
230 integer*4, dimension(nfunc) :: func_type
231 integer*4, dimension(nfunc +1) :: func_indx
232
233
234 integer*4, dimension(nhexa,9) :: con_hexa
235integer*4, dimension(max_num_ns,nl_sism) :: sour_ns
236integer*4, dimension(max_num_ne,nl_expl) :: sour_ne
237
238integer*4, dimension(:), allocatable :: node_counter_poiX,node_counter_poiY
239
240 recv_poiy, recv_poix
241
242 real*8 :: dxdx,dxdy,dxdz,dydx,dydy,dydz,dzdx,dzdy,dzdz,det_j
243 real*8 :: lambda,mu,alpha,tref,pi,lambda2, mu2, rho2
244 real*8 :: rho,ellez
245 real*8 :: l1x,l1y,l1z,l2x,l2y,l2z,area,v1,v2,v3,v4,v,term,rr
246 real*8 :: x0,y0,z0,x1,x2,x3,r1,r2,r3,f1,f2,f3,phii,theta,psi
248 real*8 :: x,y,z
249 real*8 :: csi, eta, zeta, normal_x,normal_y,normal_z
250 real*8 :: slip1_sism,slip2_sism,slip3_sism
251 real*8 :: norm1_sism,norm2_sism,norm3_sism
252 real*8 :: amp_sism
253 real*8 :: tau_sism
254 real*8 :: slip1_expl,slip2_expl,slip3_expl
255 real*8 :: norm1_expl,norm2_expl,norm3_expl
256 real*8 :: amp_expl
257 real*8 :: z_node_bottom,z_node_bottom1,z_node_bottom2
258 real*8 :: a1, a2, a3, b1, b2, b3, c1, c2, c3, w1, cost1, cost2, cost3
259 real*8 :: u1, u2, u3, u1_12, u1_13, u2_12, u2_23, u3_13, u3_23, k1
260 real*8 :: vp,vs,omega, n1, n2, n3
261
262 real*8, dimension(:), allocatable :: val_glox, val_gloy, val_gloz
263 real*8, dimension(:), allocatable :: val_locx, val_locy, val_locz
264 real*8, dimension(:), allocatable :: ct,ww
265 real*8, dimension(:), allocatable :: normal_nx_el_neun
266 real*8, dimension(:), allocatable :: normal_ny_el_neun
267 real*8, dimension(:), allocatable :: normal_nz_el_neun
268 real*8, dimension(:), allocatable :: dist_tra, dist_tra_glo
269 real*8, dimension(:), allocatable :: x_tra,y_tra,z_tra,x_tra_glo,y_tra_glo
270 real*8, dimension(:), allocatable :: xt_tra,yt_tra,zt_tra,xt_tra_glo
271 real*8, dimension(:), allocatable :: dist_tra_real
272 real*8, dimension(nnod_loc) :: xs_loc,ys_loc,zs_loc
273 real*8, dimension(nnod_loc) :: lambda_nhe, rho_nhe, mu_nhe
274 real*8, dimension(:,:,:), allocatable :: lambda_el, rho_el, mu_el,
275 real*8, dimension(nm) :: tref_mat
276 real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13, alfa21,alfa22,alfa23
277 real*8, dimension(ne_loc) :: beta11,beta12,beta13, beta21,beta22,beta23
278 real*8, dimension(ne_loc) :: gamma1,gamma2,gamma3, delta1,delta2,delta3
279 real*8, dimension(6) :: sum_facs
280
281 real*8, dimension(:,:), allocatable :: dd
282 real*8, dimension(nm,4) :: prop_mat
283 real*8, dimension(nl_dirX,4) :: val_dirx
284 real*8, dimension(nl_dirY,4) :: val_diry
285 real*8, dimension(nl_dirZ,4) :: val_dirz
286 real*8, dimension(nl_neuX,4) :: val_neux
287 real*8, dimension(nl_neuY,4) :: val_neuy
288 real*8, dimension(nl_neuZ,4) :: val_neuz
289 real*8, dimension(nl_neuN,4) :: val_neun
290 real*8, dimension(nl_poiX,4) :: val_poix
291 real*8, dimension(nl_poiY,4) :: val_poiy
292 real*8, dimension(nl_poiZ,4) :: val_poiz
293 real*8, dimension(nl_traX,4) :: val_trax
294 real*8, dimension(nl_traY,4) :: val_tray
295 real*8, dimension(nl_traZ,4) :: val_traz
296 real*8, dimension(nl_plaX,1) :: val_plax
297 real*8, dimension(nl_plaY,1) :: val_play
298 real*8, dimension(nl_plaZ,1) :: val_plaz
299 real*8, dimension(nl_sism,szsism) :: val_sism
300 real*8, dimension(nl_expl,20) :: val_expl
301 real*8, dimension(nl_forX,4) :: val_forx
302 real*8, dimension(nl_forY,4) :: val_fory
303 real*8, dimension(nl_forZ,4) :: val_forz
304 real*8, dimension(nl_forc,10) :: val_forc
305 real*8, dimension(nl_pres,10) :: val_pres
306 real*8, dimension(nl_shea,10) :: val_shea
307
308 real*8, dimension(nl_expl,6) :: facsexpl
309 real*8, dimension(nfunc,3*nnod_loc) :: fmat
310 real*8, dimension(nl_sism,6) :: facsmom
311 real*8, dimension(nl_sism,1) :: tausmom
312 real*8, dimension(3,3) :: rot
313 real*8, dimension(nfunc_data) :: func_data
314
315
316
317 pi = 4.d0*datan(1.d0);
318
319
320
321
322 length_cns = 0
323 if (nl_sism.gt.0) then
324 do isism = 1,nl_sism
325 facsmom(isism,1) = 0.0d0
326 facsmom(isism,2) = 0.0d0
327 facsmom(isism,3) = 0.0d0
328 facsmom(isism,4) = 0.0d0
329 facsmom(isism,5) = 0.0d0
330 facsmom(isism,6) = 0.0d0
331 tausmom(isism,1) = 0.0d0
332 enddo
333 endif
334
335
336
337 length_cne = 0
338 if (nl_expl.gt.0) then
339 do iexpl = 1,nl_expl
340 facsexpl(iexpl,1) = 0.0d0
341 facsexpl(iexpl,2) = 0.0d0
342 facsexpl(iexpl,3) = 0.0d0
343 facsexpl(iexpl,4) = 0.0d0
344 facsexpl(iexpl,5) = 0.0d0
345 facsexpl(iexpl,6) = 0.0d0
346 enddo
347 endif
348
349
350
351
352
353
354
355
356
357
358
359 if (nl_neun.gt.0) then
360
361
362 nelem_neun = 0
363 allocate(i4count(nnod_loc))
364 i4count = 0
365
367 nnode_neun, i4count, local_n_num)
368
369
370
371
372
373
374
375
376 do im = 1,nm
377 nn = sdeg_mat(im) +1
378
379 do ie = 1,ne_loc
380
381 if (cs_loc(cs_loc(ie -1) +0).eq.tag_mat(im)) then
382 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -
383 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1
384 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1
385 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1
386
387 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
388 nelem_neun = nelem_neun +1
389 endif
390
391 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -
392 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn
393 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn
394 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1
395
396 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
397 nelem_neun = nelem_neun +1
398 endif
399
400 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -
401 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn
402 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn
403 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1
404
405 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
406 nelem_neun = nelem_neun +1
407 endif
408
409 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn
410 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn
411 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn
412 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn
413
414 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
415 nelem_neun = nelem_neun +1
416 endif
417
418 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1
419 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn
420 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn
421 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1
422
423 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
424 nelem_neun = nelem_neun +1
425 endif
426
427 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1
428 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn
429 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn
430 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1
431
432 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
433 nelem_neun = nelem_neun +1
434 endif
435
436 endif
437 enddo
438 enddo
439
440
441
442 allocate(i4normal(nelem_neun))
443 allocate(normal_nx_el_neun(nelem_neun))
444 allocate(normal_ny_el_neun(nelem_neun))
445 allocate(normal_nz_el_neun(nelem_neun))
446
447 normal_nx_el_neun = 0.0d0
448 normal_ny_el_neun = 0.0d0
449 normal_nz_el_neun = 0.0d0
450 nelem_neun = 0
451
452 do im = 1,nm
453 nn = sdeg_mat(im) +1
454
455 do ie = 1,ne_loc
456
457 if (cs_loc(cs_loc(ie -1) +0).eq.tag_mat(im)) then
458 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1
459 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -
460 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1
461 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -
462
463 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
464
465
466
467 nelem_neun = nelem_neun +1
468
470 if(ie_surf .eq. 0) then
471 write(*,*) 'surf not found'
473 endif
474 i4normal(nelem_neun) = ie_surf
475
476 call make_normal(1,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3
477 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3
478 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3
479 normal_x,normal_y,normal_z, 0, 0 )
480
481
482 normal_nx_el_neun(nelem_neun) = normal_x
483 normal_ny_el_neun(nelem_neun) = normal_y
484 normal_nz_el_neun(nelem_neun) = normal_z
485
486
487 endif
488
489 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -
490 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn
491 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn
492 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1
493
494 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
495
496
497 nelem_neun = nelem_neun +1
498
500 if(ie_surf .eq. 0) then
501 write(*,*) 'surf not found'
503 endif
504 i4normal(nelem_neun) = ie_surf
505
506 call make_normal(2,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3)
507 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3)
508 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3)
509 normal_x,normal_y,normal_z, 0, 0 )
510
511
512 normal_nx_el_neun(nelem_neun) = normal_x
513 normal_ny_el_neun(nelem_neun) = normal_y
514 normal_nz_el_neun(nelem_neun) = normal_z
515
516 endif
517
518 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -
519 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn
520 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn
521 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1
522
523 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
524
525
526 nelem_neun = nelem_neun +1
527
529 if(ie_surf .eq. 0) then
530 write(*,*) 'surf not found'
532 endif
533 i4normal(nelem_neun) = ie_surf
534
535 call make_normal(3,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3)
536 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3)
537 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3)
538 normal_x,normal_y,normal_z, 0, 0 )
539
540
541 normal_nx_el_neun(nelem_neun) = normal_x
542 normal_ny_el_neun(nelem_neun) = normal_y
543 normal_nz_el_neun(nelem_neun) = normal_z
544
545 endif
546
547 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn
548 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn
549 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn
550 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn
551
552 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
553
554
555
556 nelem_neun = nelem_neun +1
557
559 if(ie_surf .eq. 0) then
560 write(*,*) 'surf not found'
562 endif
563
564 i4normal(nelem_neun) = ie_surf
565
566 call make_normal(4,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3)
567 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3)
568 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3)
569 normal_x,normal_y,normal_z, 0, 0 )
570
571
572 normal_nx_el_neun(nelem_neun) = normal_x
573 normal_ny_el_neun(nelem_neun) = normal_y
574 normal_nz_el_neun(nelem_neun) = normal_z
575
576 endif
577
578 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1
579 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn
580 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn
581 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1
582
583 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
584
585
586
587 nelem_neun = nelem_neun +1
588
590 if(ie_surf .eq. 0) then
591 write(*,*) 'surf not found'
593 endif
594
595
596 i4normal(nelem_neun) = ie_surf
597
598 call make_normal(5,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3)
599 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3)
600 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3)
601 normal_x,normal_y,normal_z, 0, 0 )
602
603
604 normal_nx_el_neun(nelem_neun) = normal_x
605 normal_ny_el_neun(nelem_neun) = normal_y
606 normal_nz_el_neun(nelem_neun) = normal_z
607
608 endif
609
610 ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1
611 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn
612 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn
613 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1
614
615 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4countthen
616
617
618 nelem_neun = nelem_neun +1
619
621 if(ie_surf .eq. 0) then
622 write(*,*) 'surf not found'
624 endif
625 i4normal(nelem_neun) = ie_surf
626
627 call make_normal(6,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3)
628 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3)
629 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3)
630 normal_x,normal_y,normal_z, 0, 0 )
631
632
633 normal_nx_el_neun(nelem_neun) = normal_x
634 normal_ny_el_neun(nelem_neun) = normal_y
635 normal_nz_el_neun(nelem_neun) = normal_z
636
637 endif
638
639 endif
640 enddo
641 enddo
642
643 deallocate(i4count)
644
645 endif
646
647
648
649
650
651 fmat = 0.0d0
652
653
654
655
656
657
658 if (nl_poix .gt. 0) then
659
660 allocate(val_locx(nl_poix), val_glox(nl_poix*mpi_np), ind_locx
661
662 do i = 1,nl_poix
664 val_locx(i) = dist
665 ind_locx(i) = local_n_num(node_poix(i))
666 enddo
667
668
669 call mpi_barrier(mpi_comm,mpi_ierr)
670
671 call mpi_allgather(val_locx, nl_poix, speed_double, val_glox
672 call mpi_allgather(ind_locx, nl_poix, speed_integer, ind_glox
673
674 ind_locx = 0
675
676 call get_minvalues(ind_glox, val_glox, nl_poix*mpi_np, ind_locx
677
678 do i = 1, nl_poix
679 node_poix(i) = ind_glox(ind_locx(i))
680 enddo
681
682 deallocate(ind_locx, val_locx, ind_glox, val_glox)
683 allocate(node_counter_poix(nl_poix),recv_poix(nl_poix))
684 node_counter_poix = 0;
685
686
687
688 do ie = 1,ne_loc
689 im = cs_loc(cs_loc(ie -1) +0);
690 nn = sdeg_mat(im) +1
691 do ip = 1, nl_poix
692
693 do k = 1,nn
694 do j = 1,nn
695 do i = 1,nn
696 is = nn*nn*(k -1) +nn*(j -1) +i
697 in = cs_loc(cs_loc(ie -1) +is)
698
699 if (local_n_num(in) .eq. node_poix(ip)) then
700 node_counter_poix(ip) = node_counter_poix(ip
701 endif
702 enddo
703 enddo
704 enddo
705 enddo
706
707 enddo
708
709
710
711 call mpi_barrier(mpi_comm,mpi_ierr)
712 call mpi_allreduce(node_counter_poix,recv_poix, nl_poix, speed_integer
713 node_counter_poix = recv_poix;
714 deallocate(recv_poix);
715
716 endif
717
718
719 if (nl_poiy .gt. 0) then
720
721 allocate(val_locy(nl_poiy), val_gloy(nl_poiy*mpi_np), ind_locy
722
723 do i = 1,nl_poiy
725 val_locy(i) = dist
726 ind_locy(i) = local_n_num(node_poiy(i))
727 enddo
728
729 call mpi_barrier(mpi_comm,mpi_ierr)
730
731 call mpi_allgather(val_locy, nl_poiy, speed_double, val_gloy
732 call mpi_allgather(ind_locy, nl_poiy, speed_integer, ind_gloy
733
734 call get_minvalues(ind_gloy,val_gloy, nl_poiy*mpi_np, ind_locy
735
736
737 do i = 1, nl_poiy
738 node_poiy(i) = ind_gloy(ind_locy(i))
739 enddo
740
741 deallocate(ind_locy, val_locy, ind_gloy, val_gloy)
742 allocate(node_counter_poiy(nl_poiy),recv_poiy(nl_poiy))
743 node_counter_poiy = 0;
744
745
746
747 do ie = 1,ne_loc
748 im = cs_loc(cs_loc(ie -1) +0);
749 nn = sdeg_mat(im) +1
750 do ip = 1, nl_poiy
751
752 do k = 1,nn
753 do j = 1,nn
754 do i = 1,nn
755 is = nn*nn*(k -1) +nn*(j -1) +i
756 in = cs_loc(cs_loc(ie -1) +is)
757
758 if (local_n_num(in) .eq. node_poiy(ip)) then
759 node_counter_poiy(ip) = node_counter_poiy(ip
760 endif
761 enddo
762 enddo
763 enddo
764 enddo
765
766 enddo
767
768
769
770 call mpi_barrier(mpi_comm,mpi_ierr)
771 call mpi_allreduce(node_counter_poiy,recv_poiy, nl_poiy, speed_integer
772 node_counter_poiy = recv_poiy;
773 deallocate(recv_poiy);
774
775
776 endif
777
778 if (nl_poiz .gt. 0) then
779
780 allocate(val_locz(nl_poiz), val_gloz(nl_poiz*mpi_np), ind_locz
781
782 do i = 1,nl_poiz
784 val_locz(i) = dist
785 ind_locz(i) = local_n_num(node_poiz(i))
786 enddo
787
788 call mpi_barrier(mpi_comm,mpi_ierr)
789
790 call mpi_allgather(val_locz, nl_poiz, speed_double, val_gloz
791 call mpi_allgather(ind_locz, nl_poiz, speed_integer, ind_gloz
792 call get_minvalues(ind_gloz, val_gloz, nl_poiz*mpi_np, ind_locz
793
794 do i = 1, nl_poiz
795 node_poiz(i) = ind_gloz(ind_locz(i))
796
797 enddo
798
799
800 deallocate(ind_locz, val_locz, ind_gloz, val_gloz)
801 allocate(node_counter_poiz(nl_poiz),recv_poiz(nl_poiz))
802 node_counter_poiz = 0;
803
804
805
806 do ie = 1,ne_loc
807 im = cs_loc(cs_loc(ie -1) +0);
808 nn = sdeg_mat(im) +1
809 do ip = 1, nl_poiz
810
811 do k = 1,nn
812 do j = 1,nn
813 do i = 1,nn
814 is = nn*nn*(k -1) +nn*(j -1) +i
815 in = cs_loc(cs_loc(ie -1) +is)
816
817 if (local_n_num(in) .eq. node_poiz(ip)) then
818 node_counter_poiz(ip) = node_counter_poiz(ip
819 endif
820 enddo
821 enddo
822 enddo
823 enddo
824
825 enddo
826
827
828
829 call mpi_barrier(mpi_comm,mpi_ierr)
830 call mpi_allreduce(node_counter_poiz, recv_poiz, nl_poiz, speed_integer
831 node_counter_poiz = recv_poiz;
832 deallocate(recv_poiz);
833
834 call mpi_barrier(mpi_comm,mpi_ierr)
835
836
837 endif
838
839 ne_loc = cs_loc(0) -1
840
841
842
843 do ie = 1, ne_loc
844
845 im = cs_loc(cs_loc(ie -1) +0);
846
847 nn = sdeg_mat(im) +1
848 allocate(ct(nn),ww(nn),dd(nn,nn))
850
851 rho = prop_mat(im,1)
852 lambda = prop_mat(im,2)
853 mu = prop_mat(im,3)
854
855 if (nmat_nhe.gt.0) then
856 allocate(rho_el(nn,nn,nn), lambda_el(nn,nn,nn), mu_el(nn
858 rho_nhe, lambda_nhe, mu_nhe
859 rho_el, lambda_el, mu_el
860 rho = sum(rho_el)/(nn*nn*nn)
861 lambda = sum(lambda_el)/(nn*nn*nn)
862 mu = sum(mu_el)/(nn*nn*nn)
863 deallocate(rho_el, lambda_el, mu_el, gamma_el)
864 endif
865
866 if (nl_poix.gt.0) then
867 do ip = 1,nl_poix
868 fn = 0
869 do ifun = 1,nfunc
870 if (fun_poix(ip).eq.tag_func(ifun)) fn = ifun
871 enddo
872
873 if (fn.gt.0) then
874
876
877 if (ic1 .ne. 0) then
878
879 do k = 1,nn
880 do j = 1,nn
881 do i = 1,nn
882 is = nn*nn*(k -1) +nn*(j -1) +i
883 in = cs_loc(cs_loc(ie -1) +is)
884
885 if (local_n_num(in) .eq. node_poix(ipthen
886 fmat(fn,(3*(in -1) +1)) = fmat(fn
887
888
889
890 endif
891 enddo
892 enddo
893 enddo
894 endif
895
896 endif
897 enddo
898 endif
899
900 if (nl_poiy.gt.0) then
901 do ip = 1,nl_poiy
902 fn = 0
903 do ifun = 1,nfunc
904 if (fun_poiy(ip).eq.tag_func(ifun)) fn = ifun
905 enddo
906
907 if (fn.gt.0) then
908
910 if (ic1 .ne. 0) then
911 do k = 1,nn
912 do j = 1,nn
913 do i = 1,nn
914
915 is = nn*nn*(k -1) +nn*(j -1) +i
916 in = cs_loc(cs_loc(ie -1) +is)
917
918 if (local_n_num(in) .eq. node_poiy(ipthen
919 fmat(fn,(3*(in -1) +2)) = fmat(fn
920 endif
921 enddo
922 enddo
923 enddo
924 endif
925 endif
926 enddo
927 endif
928
929 if (nl_poiz.gt.0) then
930 do ip = 1,nl_poiz
931 fn = 0
932 do ifun = 1,nfunc
933 if (fun_poiz(ip).eq.tag_func(ifun)) fn = ifun
934 enddo
935
936 if (fn.gt.0) then
937
939 if (ic1 .ne. 0) then
940
941 do k = 1,nn
942 do j = 1,nn
943 do i = 1,nn
944
945 is = nn*nn*(k -1) +nn*(j -1) +i
946 in = cs_loc(cs_loc(ie -1) +is)
947
948 if (local_n_num(in) .eq. node_poiz(ip)) then
949
950 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in
951 + val_poiz(ip,4)/node_counter_poiz
952
953 endif
954
955 enddo
956 enddo
957 enddo
958 endif
959 endif
960 enddo
961 endif
962
963
964
965
966
967
968 if (nl_plax.gt.0) then
969
970 do ipl = 1,nl_plax
971
972 if (tag_mat(im).eq.tag_plax(ipl)) then
973
974
975 sum_node_bottom1 = 0
976 sum_node_bottom2 = 0
977 last_node_bottom1 = 0
978 last_node_bottom2 = 0
979
980 do i = 1,8
981 ie_glob = local_el_num(ie)
982
985
986 if(ic1 .ne. 0 .and. ic2 .ne. 0) then
987 if (dabs(zs_loc(ic1) - zs_loc(ic2)) .le. 1.d-4) then
988 last_node_bottom1 = i
989 sum_node_bottom1 = sum_node_bottom1 + i
990 z_node_bottom1 = zs_loc(ic1)
991
992 else
993 last_node_bottom2 = i
994 sum_node_bottom2 = sum_node_bottom2 + i
995 z_node_bottom2 = zs_loc(ic2)
996
997 endif
998 endif
999
1000 enddo
1001
1002 if (z_node_bottom1.lt.z_node_bottom2) then
1003 sum_node_bottom = sum_node_bottom1
1004 sum_node_bottom_first3 = sum_node_bottom1 - last_node_bottom1
1005 else
1006 sum_node_bottom = sum_node_bottom2
1007 sum_node_bottom_first3 = sum_node_bottom2 - last_node_bottom2
1008 endif
1009
1010
1011 c=dsqrt(mu*rho)
1012 ellez=dabs(z_node_bottom1 - z_node_bottom2)
1013
1014 fn = 0
1015 do ifun = 1,nfunc
1016 if (fun_plax(ipl).eq.tag_func(ifun)) fn = ifun
1017 enddo
1018 if (fn.gt.0) then
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079 if ((sum_node_bottom.eq.10).or.(sum_node_bottom.eq.then
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148 if (sum_node_bottom.eq.10) then
1149
1150 sit = 1
1151
1152
1153
1154 else
1155
1156 sit = 2
1157
1158 endif
1159
1160
1161 if (mod(nn-1,2).eq.0) then
1162 k = ((nn-1)/2)+1
1163
1164
1165 else
1166 k = (int((nn-1)/2))+num_sit(sit)
1167 endif
1168
1169
1170 do i = 1,nn
1171 do j = 1,nn
1172 dxdx = alfa11(ie) +beta12(ie)*ct(k
1173 + beta13(ie)*ct(j) +gamma1
1174 dydx = alfa21(ie) +beta22(ie)*ct(k
1175 + beta23(ie)*ct(j) +gamma2
1176 dzdx = alfa31(ie) +beta32(ie)*ct(k
1177 + beta33(ie)*ct(j) +gamma3
1178 dxdy = alfa12(ie) +beta11(ie)*ct(k
1179 + beta13(ie)*ct(i) +gamma1
1180 dydy = alfa22(ie) +beta21(ie)*ct(k
1181 + beta23(ie)*ct(i) +gamma2
1182 dzdy = alfa32(ie) +beta31(ie)*ct(k
1183 + beta33(ie)*ct(i) +gamma3
1184 dxdz = alfa13(ie) +beta11(ie)*ct(j
1185 + beta12(ie)*ct(i) +gamma1
1186 dydz = alfa23(ie) +beta21(ie)*ct(j
1187 + beta22(ie)*ct(i) +gamma2
1188 dzdz = alfa33(ie) +beta31(ie)*ct(j
1189 + beta32(ie)*ct(i) +gamma3
1190 det_j = dxdz * (dydx*dzdy - dzdx*dydy
1191 - dydz * (dxdx*dzdy - dzdx*dxdy
1192 + dzdz * (dxdx*dydy - dydx*dxdy
1193 is = nn*nn*(k -1) +nn*(j -1) +i
1194 in = cs_loc(cs_loc(ie -1) + is)
1195
1196 term = 4 * c * val_plax(ipl,1) * ww
1197 fmat(fn,(3*(in -1) +1)) = fmat(fn
1198
1199 prova = 3*(in -1) +1
1200
1201 enddo
1202 enddo
1203
1204
1205
1206
1207
1208
1209 elseif ((sum_node_bottom.eq.14).or.(sum_node_bottom.eq.then
1210
1211
1212
1213
1214 if (sum_node_bottom.eq.14) then
1215
1216 sit = 1
1217
1218
1219
1220 else
1221
1222 sit = 2
1223
1224 endif
1225
1226
1227 if (mod(nn-1,2).eq.0) then
1228 j = ((nn-1)/2)+1
1229
1230
1231 else
1232 j = (int((nn-1)/2))+num_sit(sit)
1233 endif
1234
1235 do i = 1,nn
1236 do k = 1,nn
1237 dxdx = alfa11(ie) +beta12(ie)*ct(k)
1238 + beta13(ie)*ct(j) +gamma1(ie
1239 dydx = alfa21(ie) +beta22(ie)*ct(k)
1240 + beta23(ie)*ct(j) +gamma2(ie
1241 dzdx = alfa31(ie) +beta32(ie)*ct(k)
1242 + beta33(ie)*ct(j) +gamma3(ie
1243 dxdy = alfa12(ie) +beta11(ie)*ct(k)
1244 + beta13(ie)*ct(i) +gamma1(ie
1245 dydy = alfa22(ie) +beta21(ie)*ct(k)
1246 + beta23(ie)*ct(i) +gamma2(ie
1247 dzdy = alfa32(ie) +beta31(ie)*ct(k)
1248 + beta33(ie)*ct(i) +gamma3(ie
1249 dxdz = alfa13(ie) +beta11(ie)*ct(j)
1250 + beta12(ie)*ct(i) +gamma1(ie
1251 dydz = alfa23(ie) +beta21(ie)*ct(j)
1252 + beta22(ie)*ct(i) +gamma2(ie
1253 dzdz = alfa33(ie) +beta31(ie)*ct(j)
1254 + beta32(ie)*ct(i) +gamma3(ie
1255 det_j = dxdz * (dydx*dzdy - dzdx*dydy
1256 - dydz * (dxdx*dzdy - dzdx*dxdy
1257 + dzdz * (dxdx*dydy - dydx*dxdy
1258
1259 is = nn*nn*(k -1) +nn*(j -1) +i
1260 in = cs_loc(cs_loc(ie -1) +is)
1261
1262 term = 4 * c * val_plax(ipl,1) * ww
1263
1264 prova = 3*(in -1) +1
1265
1266
1267 fmat(fn,(3*(in -1) +1)) = fmat(fn,(
1268 enddo
1269 enddo
1270
1271
1272
1273
1274
1275 elseif (sum_node_bottom.eq.18) then
1276
1277
1278
1279 if ((sum_node_bottom_first3.eq.10) .or.(sum_node_bottom_first3.eq.
1280 .or.(sum_node_bottom_first3.eq.14)) then
1281
1282 sit = 1
1283
1284
1285
1286 else
1287
1288 sit = 2
1289
1290 endif
1291
1292
1293 if (mod(nn-1,2).eq.0) then
1294 i = ((nn-1)/2)+1
1295
1296
1297 else
1298 i = (int((nn-1)/2))+num_sit(sit)
1299 endif
1300
1301 do j = 1,nn
1302 do k = 1,nn
1303 dxdx = alfa11(ie) +beta12(ie)*ct(k
1304 + beta13(ie)*ct(j) +gamma1(ie
1305 dydx = alfa21(ie) +beta22(ie)*ct(k
1306 + beta23(ie)*ct(j) +gamma2(ie
1307 dzdx = alfa31(ie) +beta32(ie)*ct(k
1308 + beta33(ie)*ct(j) +gamma3(ie
1309 dxdy = alfa12(ie) +beta11(ie)*ct(k
1310 + beta13(ie)*ct(i) +gamma1(ie
1311 dydy = alfa22(ie) +beta21(ie)*ct(k
1312 + beta23(ie)*ct(i) +gamma2(ie
1313 dzdy = alfa32(ie) +beta31(ie)*ct(k
1314 + beta33(ie)*ct(i) +gamma3(ie
1315 dxdz = alfa13(ie) +beta11(ie)*ct(j
1316 + beta12(ie)*ct(i) +gamma1(ie
1317 dydz = alfa23(ie) +beta21(ie)*ct(j
1318 + beta22(ie)*ct(i) +gamma2(ie
1319 dzdz = alfa33(ie) +beta31(ie)*ct(j
1320 + beta32(ie)*ct(i) +gamma3(ie
1321 det_j = dxdz * (dydx*dzdy - dzdx*dydy
1322 - dydz * (dxdx*dzdy
1323 + dzdz * (dxdx*dydy
1324
1325 is = nn*nn*(k -1) +nn*(j -1) +i
1326 in = cs_loc(cs_loc(ie -1) +is)
1327
1328 term = 4 * c * val_plax(ipl,1) * ww
1329 fmat(fn,(3*(in -1) +1)) = fmat(fn,
1330
1331
1332 enddo
1333 enddo
1334
1335 endif
1336
1337
1338 endif
1339 endif
1340
1341 enddo
1342 endif
1343
1344 if (nl_play.gt.0) then
1345 do ipl = 1,nl_play
1346
1347 if (tag_mat(im).eq.tag_play(ipl)) then
1348
1349
1350 sum_node_bottom1 = 0
1351 sum_node_bottom2 = 0
1352 last_node_bottom1 = 0
1353 last_node_bottom2 = 0
1354
1355 do i = 1,8
1356 ie_glob = local_el_num(ie)
1359
1360 if(ic1 .ne. 0 .and. ic2 .ne. 0) then
1361 if (dabs(zs_loc(ic1) - zs_loc(ic2)) .le.then
1362
1363 last_node_bottom1 = i
1364 sum_node_bottom1 = sum_node_bottom1
1365 z_node_bottom1 = zs_loc(ic1)
1366 else
1367 last_node_bottom2 = i
1368 sum_node_bottom2 = sum_node_bottom2
1369 z_node_bottom2 = zs_loc(ic2)
1370 endif
1371 endif
1372 enddo
1373
1374 if (z_node_bottom1.lt.z_node_bottom2) then
1375 sum_node_bottom = sum_node_bottom1
1376 sum_node_bottom_first3 = sum_node_bottom1
1377 else
1378 sum_node_bottom = sum_node_bottom2
1379 sum_node_bottom_first3 = sum_node_bottom2
1380 endif
1381
1382
1383
1384
1385
1386 c=dsqrt(mu*rho)
1387 ellez=dabs(z_node_bottom1 - z_node_bottom2)
1388
1389 fn = 0
1390 do ifun = 1,nfunc
1391 if (fun_play(ipl).eq.tag_func(ifun)) fn = ifun
1392 enddo
1393 if (fn.gt.0) then
1394
1395
1396
1397 if ((sum_node_bottom.eq.10).or.(sum_node_bottom.eq.then
1398
1399
1400
1401
1402 if (sum_node_bottom.eq.10) then
1403
1404 sit = 1
1405
1406
1407
1408 else
1409
1410 sit = 2
1411
1412 endif
1413
1414
1415 if (mod(nn-1,2).eq.0) then
1416 k = ((nn-1)/2)+1
1417
1418
1419 else
1420 k = (int((nn-1)/2))+num_sit(sit)
1421 endif
1422
1423 do i = 1,nn
1424 do j = 1,nn
1425 dxdx = alfa11(ie) +beta12(ie)*ct(k
1426 + beta13(ie)*ct(j) +gamma1(ie
1427 dydx = alfa21(ie) +beta22(ie)*ct(k
1428 + beta23(ie)*ct(j) +gamma2(ie
1429 dzdx = alfa31(ie) +beta32(ie)*ct(k
1430 + beta33(ie)*ct(j) +gamma3(ie
1431 dxdy = alfa12(ie) +beta11(ie)*ct(k
1432 + beta13(ie)*ct(i) +gamma1(ie
1433 dydy = alfa22(ie) +beta21(ie)*ct(k
1434 + beta23(ie)*ct(i) +gamma2(ie
1435 dzdy = alfa32(ie) +beta31(ie)*ct(k
1436 + beta33(ie)*ct(i) +gamma3(ie
1437 dxdz = alfa13(ie) +beta11(ie)*ct(j
1438 + beta12(ie)*ct(i) +gamma1(ie
1439 dydz = alfa23(ie) +beta21(ie)*ct(j
1440 + beta22(ie)*ct(i) +gamma2(ie
1441 dzdz = alfa33(ie) +beta31(ie)*ct(j
1442 + beta32(ie)*ct(i) +gamma3(ie
1443
1444 det_j = dxdz * (dydx*dzdy - dzdx*dydy
1445 - dydz * (dxdx*dzdy - dzdx*dxdy
1446 + dzdz * (dxdx*dydy - dydx*dxdy
1447
1448 is = nn*nn*(k -1) +nn*(j -1) +i
1449 in = cs_loc(cs_loc(ie -1) +is)
1450
1451 term = 4 * c * val_play(ipl,1) * ww
1452 fmat(fn,(3*(in -1) +2)) = fmat(fn,
1453 enddo
1454 enddo
1455
1456
1457
1458
1459
1460
1461 elseif ((sum_node_bottom.eq.14).or.(sum_node_bottom.eq.then
1462
1463
1464
1465 if (sum_node_bottom.eq.14) then
1466
1467 sit = 1
1468
1469
1470
1471 else
1472
1473 sit = 2
1474
1475 endif
1476
1477
1478 if (mod(nn-1,2).eq.0) then
1479 j = ((nn-1)/2)+1
1480
1481
1482 else
1483 j = (int((nn-1)/2))+num_sit(sit)
1484 endif
1485
1486 do i = 1,nn
1487 do k = 1,nn
1488
1489 dxdx = alfa11(ie) +beta12(ie)*ct(k)
1490 + beta13(ie)*ct(j) +gamma1(ie)*ct
1491 dydx = alfa21(ie) +beta22(ie)*ct(k)
1492 + beta23(ie)*ct(j) +gamma2(ie)*ct
1493 dzdx = alfa31(ie) +beta32(ie)*ct(k)
1494 + beta33(ie)*ct(j) +gamma3(ie)*ct
1495 dxdy = alfa12(ie) +beta11(ie)*ct(k)
1496 + beta13(ie)*ct(i) +gamma1(ie)*ct
1497 dydy = alfa22(ie) +beta21(ie)*ct(k)
1498 + beta23(ie)*ct(i) +gamma2(ie)*ct
1499 dzdy = alfa32(ie) +beta31(ie)*ct(k)
1500 + beta33(ie)*ct(i) +gamma3(ie)*ct
1501 dxdz = alfa13(ie) +beta11(ie)*ct(j)
1502 + beta12(ie)*ct(i) +gamma1(ie)*ct
1503 dydz = alfa23(ie) +beta21(ie)*ct(j)
1504 + beta22(ie)*ct(i) +gamma2(ie)*ct
1505 dzdz = alfa33(ie) +beta31(ie)*ct(j)
1506 + beta32(ie)*ct(i) +gamma3(ie)*ct
1507 det_j = dxdz * (dydx*dzdy - dzdx*dydy
1508 - dydz * (dxdx*dzdy - dzdx*dxdy
1509 + dzdz * (dxdx*dydy - dydx*dxdy
1510
1511 is = nn*nn*(k -1) +nn*(j -1) +i
1512 in = cs_loc(cs_loc(ie -1) +is)
1513
1514 term = 4 * c * val_play(ipl,1) * ww(i
1515 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3
1516 enddo
1517 enddo
1518
1519
1520
1521
1522
1523 elseif (sum_node_bottom.eq.18) then
1524
1525
1526
1527 if ((sum_node_bottom_first3.eq.10) &
1528 .or.(sum_node_bottom_first3.eq.13) &
1529 .or.(sum_node_bottom_first3.eq.14)) then
1530
1531 sit = 1
1532
1533
1534
1535 else
1536
1537 sit = 2
1538
1539 endif
1540
1541
1542 if (mod(nn-1,2).eq.0) then
1543 i = ((nn-1)/2)+1
1544
1545
1546 else
1547 i = (int((nn-1)/2))+num_sit(sit)
1548 endif
1549
1550 do j = 1,nn
1551 do k = 1,nn
1552
1553 dxdx = alfa11(ie) +beta12(ie)*ct(k)
1554 + beta13(ie)*ct(j) +gamma1(ie)*ct
1555 dydx = alfa21(ie) +beta22(ie)*ct(k)
1556 + beta23(ie)*ct(j) +gamma2(ie)*ct
1557 dzdx = alfa31(ie) +beta32(ie)*ct(k)
1558 + beta33(ie)*ct(j) +gamma3(ie)*ct
1559 dxdy = alfa12(ie) +beta11(ie)*ct(k)
1560 + beta13(ie)*ct(i) +gamma1(ie)*ct
1561 dydy = alfa22(ie) +beta21(ie)*ct(k)
1562 + beta23(ie)*ct(i) +gamma2(ie)*ct
1563 dzdy = alfa32(ie) +beta31(ie)*ct(k)
1564 + beta33(ie)*ct(i) +gamma3(ie)*ct
1565 dxdz = alfa13(ie) +beta11(ie)*ct(j)
1566 + beta12(ie)*ct(i) +gamma1(ie)*ct
1567 dydz = alfa23(ie) +beta21(ie)*ct(j)
1568 + beta22(ie)*ct(i) +gamma2(ie)*ct
1569 dzdz = alfa33(ie) +beta31(ie)*ct(j)
1570 + beta32(ie)*ct(i) +gamma3(ie)*ct
1571 det_j = dxdz * (dydx*dzdy - dzdx*dydy
1572 - dydz * (dxdx*dzdy - dzdx*dxdy
1573 + dzdz * (dxdx*dydy - dydx*dxdy
1574
1575 is = nn*nn*(k -1) +nn*(j -1) +i
1576 in = cs_loc(cs_loc(ie -1) +is)
1577
1578 term = 4 * c * val_play(ipl,1) * ww(j
1579 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3
1580
1581 enddo
1582 enddo
1583
1584 endif
1585
1586
1587 endif
1588 endif
1589
1590 enddo
1591 endif
1592
1593 if (nl_plaz.gt.0) then
1594 do ipl = 1,nl_plaz
1595
1596 if (tag_mat(im).eq.tag_plaz(ipl)) then
1597
1598
1599 sum_node_bottom1 = 0
1600 sum_node_bottom2 = 0
1601 last_node_bottom1 = 0
1602 last_node_bottom2 = 0
1603
1604 do i = 1,8
1605 ie_glob = local_el_num(ie)
1608
1609 if(ic1 .ne. 0 .and. ic2 .ne. 0) then
1610 if (dabs(zs_loc(ic1) - zs_loc(ic2)) .le.then
1611 last_node_bottom1 = i
1612 sum_node_bottom1 = sum_node_bottom1
1613 z_node_bottom1 = zs_loc(ic1)
1614 else
1615 last_node_bottom2 = i
1616 sum_node_bottom2 = sum_node_bottom2
1617 z_node_bottom2 = zs_loc(ic2)
1618 endif
1619 endif
1620 enddo
1621
1622 if (z_node_bottom1.lt.z_node_bottom2) then
1623 sum_node_bottom = sum_node_bottom1
1624 sum_node_bottom_first3 = sum_node_bottom1
1625 else
1626 sum_node_bottom = sum_node_bottom2
1627 sum_node_bottom_first3 = sum_node_bottom2
1628 endif
1629
1630
1631
1632
1633
1634 c=dsqrt((lambda+2*mu)*rho)
1635 ellez=dabs(z_node_bottom1 - z_node_bottom2)
1636
1637 fn = 0
1638 do ifun = 1,nfunc
1639 if (fun_plaz(ipl).eq.tag_func(ifun)) fn = ifun
1640 enddo
1641 if (fn.gt.0) then
1642
1643
1644
1645 if ((sum_node_bottom.eq.10).or.(sum_node_bottom.eq.then
1646
1647
1648
1649
1650 if (sum_node_bottom.eq.10) then
1651
1652 sit = 1
1653
1654
1655
1656 else
1657
1658 sit = 2
1659
1660 endif
1661
1662
1663 if (mod(nn-1,2).eq.0) then
1664 k = ((nn-1)/2)+1
1665
1666
1667 else
1668 k = (int((nn-1)/2))+num_sit(sit)
1669 endif
1670
1671 do i = 1,nn
1672 do j = 1,nn
1673 dxdx = alfa11(ie) +beta12(ie)*ct
1674 + beta13(ie)*ct(j) +gamma1(ie
1675 dydx = alfa21(ie) +beta22(ie)*ct
1676 + beta23(ie)*ct(j) +gamma2(ie
1677 dzdx = alfa31(ie) +beta32(ie)*ct
1678 + beta33(ie)*ct(j) +gamma3(ie
1679 dxdy = alfa12(ie) +beta11(ie)*ct
1680 + beta13(ie)*ct(i) +gamma1(ie
1681 dydy = alfa22(ie) +beta21(ie)*ct
1682 + beta23(ie)*ct(i) +gamma2(ie
1683 dzdy = alfa32(ie) +beta31(ie)*ct
1684 + beta33(ie)*ct(i) +gamma3(ie
1685 dxdz = alfa13(ie) +beta11(ie)*ct
1686 + beta12(ie)*ct(i) +gamma1(ie
1687 dydz = alfa23(ie) +beta21(ie)*ct
1688 + beta22(ie)*ct(i) +gamma2(ie
1689 dzdz = alfa33(ie) +beta31(ie)*ct
1690 + beta32(ie)*ct(i) +gamma3(ie
1691 det_j = dxdz * (dydx*dzdy - dzdx
1692 - dydz * (dxdx*dzdy - dzdx
1693 + dzdz * (dxdx*dydy - dydx
1694
1695 is = nn*nn*(k -1) +nn*(j -1) +i
1696 in = cs_loc(cs_loc(ie -1) +is)
1697
1698 term = 4 * c * val_plaz(ipl,1)
1699 fmat(fn,(3*(in -1) +3)) = fmat(fn
1700 enddo
1701 enddo
1702
1703
1704
1705
1706
1707
1708 elseif ((sum_node_bottom.eq.14).or.(sum_node_bottom.eq.then
1709
1710
1711
1712 if (sum_node_bottom.eq.14) then
1713
1714 sit = 1
1715
1716
1717
1718 else
1719
1720 sit = 2
1721
1722 endif
1723
1724
1725 if (mod(nn-1,2).eq.0) then
1726 j = ((nn-1)/2)+1
1727
1728
1729 else
1730 j = (int((nn-1)/2))+num_sit(sit)
1731 endif
1732
1733 do i = 1,nn
1734 do k = 1,nn
1735 dxdx = alfa11(ie) +beta12(ie)*ct
1736 + beta13(ie)*ct(j) +gamma1(ie
1737 dydx = alfa21(ie) +beta22(ie)*ct
1738 + beta23(ie)*ct(j) +gamma2(ie
1739 dzdx = alfa31(ie) +beta32(ie)*ct
1740 + beta33(ie)*ct(j) +gamma3(ie
1741
1742 dxdy = alfa12(ie) +beta11(ie)*ct
1743 + beta13(ie)*ct(i) +gamma1(ie
1744 dydy = alfa22(ie) +beta21(ie)*ct
1745 + beta23(ie)*ct(i) +gamma2(ie
1746 dzdy = alfa32(ie) +beta31(ie)*ct
1747 + beta33(ie)*ct(i) +gamma3(ie
1748
1749 dxdz = alfa13(ie) +beta11(ie)*ct
1750 + beta12(ie)*ct(i) +gamma1(ie
1751 dydz = alfa23(ie) +beta21(ie)*ct
1752 + beta22(ie)*ct(i) +gamma2(ie
1753 dzdz = alfa33(ie) +beta31(ie)*ct
1754 + beta32(ie)*ct(i) +gamma3(ie
1755
1756 det_j = dxdz * (dydx*dzdy - dzdx
1757 - dydz * (dxdx*dzdy - dzdx
1758 + dzdz * (dxdx*dydy - dydx
1759
1760 is = nn*nn*(k -1) +nn*(j -1) +i
1761 in = cs_loc(cs_loc(ie -1) +is)
1762
1763 term = 4 * c * val_plaz(ipl,1)
1764 fmat(fn,(3*(in -1) +3)) = fmat(fn
1765 enddo
1766 enddo
1767
1768
1769
1770
1771
1772 elseif (sum_node_bottom.eq.18) then
1773
1774
1775
1776 if ((sum_node_bottom_first3.eq.10) &
1777 .or.(sum_node_bottom_first3.eq.13) &
1778 .or.(sum_node_bottom_first3.eq.14)) then
1779
1780 sit = 1
1781
1782
1783
1784 else
1785
1786 sit = 2
1787
1788 endif
1789
1790
1791 if (mod(nn-1,2).eq.0) then
1792 i = ((nn-1)/2)+1
1793
1794
1795 else
1796 i = (int((nn-1)/2))+num_sit(sit)
1797 endif
1798
1799 do j = 1,nn
1800 do k = 1,nn
1801
1802 dxdx = alfa11(ie) +beta12(ie)*ct(k)
1803 + beta13(ie)*ct(j) +gamma1(ie)
1804 dydx = alfa21(ie) +beta22(ie)*ct(k)
1805 + beta23(ie)*ct(j) +gamma2(ie)
1806 dzdx = alfa31(ie) +beta32(ie)*ct(k)
1807 + beta33(ie)*ct(j) +gamma3(ie)
1808
1809 dxdy = alfa12(ie) +beta11(ie)*ct(k)
1810 + beta13(ie)*ct(i) +gamma1(ie)
1811 dydy = alfa22(ie) +beta21(ie)*ct(k)
1812 + beta23(ie)*ct(i) +gamma2(ie)
1813 dzdy = alfa32(ie) +beta31(ie)*ct(k)
1814 + beta33(ie)*ct(i) +gamma3(ie)
1815
1816 dxdz = alfa13(ie) +beta11(ie)*ct(j)
1817 + beta12(ie)*ct(i) +gamma1(ie)
1818 dydz = alfa23(ie) +beta21(ie)*ct(j)
1819 + beta22(ie)*ct(i) +gamma2(ie)
1820 dzdz = alfa33(ie) +beta31(ie)*ct(j)
1821 + beta32(ie)*ct(i) +gamma3(ie)
1822
1823 det_j = dxdz * (dydx*dzdy - dzdx*dydy
1824 - dydz * (dxdx*dzdy - dzdx*dxdy
1825 + dzdz * (dxdx*dydy - dydx*dxdy
1826
1827 is = nn*nn*(k -1) +nn*(j -1) +i
1828 in = cs_loc(cs_loc(ie -1) +is)
1829
1830 term = 4 * c * val_plaz(ipl,1)
1831 fmat(fn,(3*(in -1) +3)) = fmat(fn
1832 enddo
1833 enddo
1834
1835 endif
1836
1837
1838 endif
1839 endif
1840
1841 enddo
1842 endif
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855 if (nl_sism.gt.0) then
1856
1857 do isism = 1, nl_sism
1858 do ipsism = 1, num_ns(isism)
1859
1860 fn = 0
1861
1862 do ifun = 1,nfunc
1863 if (fun_sism(isism).eq.tag_func(ifun)) fn = ifun
1864 enddo
1865
1866 if (fn.gt.0) then
1867 do k = 1,nn
1868 do j = 1,nn
1869 do i = 1,nn
1870 dxdx = alfa11(ie) +beta12(ie)*ct(k) &
1871 + beta13(ie)*ct(j) +gamma1(ie)*ct
1872 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1873 + beta23(ie)*ct(j) +gamma2(ie)*ct
1874 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1875 + beta33(ie)*ct(j) +gamma3(ie)*ct
1876
1877 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1878 + beta13(ie)*ct(i) +gamma1(ie)*ct
1879 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1880 + beta23(ie)*ct(i) +gamma2(ie)*ct
1881 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1882 + beta33(ie)*ct(i) +gamma3(ie)*ct
1883
1884 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1885 + beta12(ie)*ct(i) +gamma1(ie)*ct
1886 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1887 + beta22(ie)*ct(i) +gamma2(ie)*ct
1888 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1889 + beta32(ie)*ct(i) +gamma3(ie)*ct
1890
1891 det_j = dxdz * (dydx*dzdy - dzdx*dydy)
1892 - dydz * (dxdx*dzdy - dzdx*dxdy)
1893 + dzdz * (dxdx*dydy - dydx*dxdy)
1894
1895 is = nn*nn*(k -1) +nn*(j -1) +i
1896 in = cs_loc(cs_loc(ie -1) +is)
1897
1898
1899 if (local_n_num(in) .eq. sour_ns(ipsismthen
1900 facsmom(isism,1) = facsmom(isism
1901 facsmom(isism,2) = facsmom(isism
1902 facsmom(isism,3) = facsmom(isism
1903 facsmom(isism,4) = facsmom(isism
1904 facsmom(isism,5) = facsmom(isism
1905 facsmom(isism,6) = facsmom(isism
1906
1907 length_cns = length_cns + 1
1908
1909 endif
1910 enddo
1911 enddo
1912 enddo
1913 endif
1914
1915 enddo
1916
1917 enddo
1918
1919 endif
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931 if (nl_expl.gt.0) then
1932 do iexpl = 1,nl_expl
1933 do ipexpl = 1,num_ne(iexpl)
1934
1935 fn = 0
1936 do ifun = 1,nfunc
1937 if (fun_expl(iexpl).eq.tag_func(ifun)) fn = ifun
1938 enddo
1939
1940 if (fn.gt.0) then
1941 do k = 1,nn
1942 do j = 1,nn
1943 do i = 1,nn
1944 dxdx = alfa11(ie) +beta12(ie)*ct(k) &
1945 + beta13(ie)*ct(j) +gamma1(ie)*ct
1946 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1947 + beta23(ie)*ct(j) +gamma2(ie)*ct
1948 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1949 + beta33(ie)*ct(j) +gamma3(ie)*ct
1950
1951 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1952 + beta13(ie)*ct(i) +gamma1(ie)*ct
1953 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1954 + beta23(ie)*ct(i) +gamma2(ie)*ct
1955 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1956 + beta33(ie)*ct(i) +gamma3(ie)*ct
1957
1958 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1959 + beta12(ie)*ct(i) +gamma1(ie)*ct
1960 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1961 + beta22(ie)*ct(i) +gamma2(ie)*ct
1962 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1963 + beta32(ie)*ct(i) +gamma3(ie)*ct
1964
1965 det_j = dxdz * (dydx*dzdy - dzdx*dydy)
1966 - dydz * (dxdx*dzdy - dzdx*dxdy)
1967 + dzdz * (dxdx*dydy - dydx*dxdy)
1968
1969 is = nn*nn*(k -1) +nn*(j -1) +i
1970 in = cs_loc(cs_loc(ie -1) +is)
1971
1972 if (local_n_num(in) .eq. sour_ne(ipexplthen
1973 facsexpl(iexpl,1) = facsexpl(iexpl
1974 facsexpl(iexpl,2) = facsexpl(iexpl
1975 facsexpl(iexpl,3) = facsexpl(iexpl
1976 facsexpl(iexpl,4) = facsexpl(iexpl
1977 facsexpl(iexpl,5) = facsexpl(iexpl
1978 facsexpl(iexpl,6) = facsexpl(iexpl
1979
1980 length_cne = length_cne + 1
1981
1982 endif
1983 enddo
1984 enddo
1985 enddo
1986 endif
1987
1988 enddo
1989
1990 enddo
1991
1992
1993 endif
1994
1995
1996
1997
1998
1999
2000 deallocate(ct,dd,ww)
2001 enddo
2002
2003
2004
2005
2006
2007
2008
2009 if (nl_sism.gt.0) then
2010 do isism = 1,nl_sism
2011 call mpi_barrier(mpi_comm,mpi_ierr)
2012 call mpi_allreduce(facsmom(isism,:),sum_facs, 6, speed_double
2013 facsmom(isism,:) = sum_facs
2014 enddo
2015
2016 endif
2017
2018 if (nl_sism.gt.0) then
2019 if (srcmodflag.eq.0) then
2020 val_st = 12
2021 elseif (srcmodflag.eq.1) then
2022 val_st = 6
2023 endif
2024 do isism = 1,nl_sism
2025 slip1_sism = val_sism(isism,val_st + 1)
2026 slip2_sism = val_sism(isism,val_st + 2)
2027 slip3_sism = val_sism(isism,val_st + 3)
2028
2029 norm1_sism = val_sism(isism,val_st + 4)
2030 norm2_sism = val_sism(isism,val_st + 5)
2031 norm3_sism = val_sism(isism,val_st + 6)
2032
2033 amp_sism = val_sism(isism,val_st + 8)
2034
2035 tau_sism = val_sism(isism,val_st + 9)
2036
2037
2038
2039
2040 if(facsmom(isism,1) .ne. 0) &
2041 facsmom(isism,1) = 1/facsmom(isism
2042 * (slip1_sism*norm1_sism+slip1_sism
2043 * amp_sism
2044 if(facsmom(isism,2) .ne. 0) &
2045 facsmom(isism,2) = 1/facsmom(isism
2046 * (slip2_sism*norm2_sism+slip2_sism
2047 * amp_sism
2048 if(facsmom(isism,3) .ne. 0) &
2049 facsmom(isism,3) = 1/facsmom(isism
2050 * (slip3_sism*norm3_sism+slip3_sism
2051 * amp_sism
2052 if(facsmom(isism,4) .ne. 0) &
2053 facsmom(isism,4) = 1/facsmom(isism
2054 * (slip2_sism*norm3_sism+slip3_sism
2055 * amp_sism
2056 if(facsmom(isism,5) .ne. 0) &
2057 facsmom(isism,5) = 1/facsmom(isism
2058 * (slip1_sism*norm3_sism+slip3_sism
2059 * amp_sism
2060 if(facsmom(isism,6) .ne. 0) &
2061 facsmom(isism,6) = 1/facsmom(isism
2062 * (slip1_sism*norm2_sism+slip2_sism
2063 * amp_sism
2064
2065 tausmom(isism,1) = tau_sism
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094 enddo
2095
2096 endif
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108 if (nl_expl.gt.0) then
2109 do iexpl = 1,nl_expl
2110
2111 call mpi_barrier(mpi_comm,mpi_ierr)
2112 call mpi_allreduce(facsexpl(iexpl,:),sum_facs, 6, speed_double
2113 facsexpl(iexpl,:) = sum_facs
2114 enddo
2115 endif
2116
2117
2118
2119 if (nl_expl.gt.0) then
2120
2121 do iexpl = 1,nl_expl
2122 slip1_expl = val_expl(iexpl,13)
2123 slip2_expl = val_expl(iexpl,14)
2124 slip3_expl = val_expl(iexpl,15)
2125
2126 norm1_expl = val_expl(iexpl,16)
2127 norm2_expl = val_expl(iexpl,17)
2128 norm3_expl = val_expl(iexpl,18)
2129
2130 amp_expl = val_expl(iexpl,20)
2131
2132 facsexpl(iexpl,1) = 1/facsexpl(iexpl,1) &
2133 * slip1_expl &
2134 * amp_expl
2135 facsexpl(iexpl,2) = 1/facsexpl(iexpl,2) &
2136 * slip2_expl &
2137 * amp_expl
2138 facsexpl(iexpl,3) = 1/facsexpl(iexpl,3) &
2139 * slip3_expl &
2140 * amp_expl
2141 facsexpl(iexpl,4) = 1/facsexpl(iexpl,4) &
2142 * norm1_expl &
2143 * amp_expl
2144 facsexpl(iexpl,5) = 1/facsexpl(iexpl,5) &
2145 * norm2_expl &
2146 * amp_expl
2147 facsexpl(iexpl,6) = 1/facsexpl(iexpl,6) &
2148 * norm3_expl &
2149 * amp_expl
2150 enddo
2151 endif
2152
2153
2154
2155
2156
2157 ne_loc = cs_loc(0) -1
2158
2159 do ie = 1,ne_loc
2160 im = cs_loc(cs_loc(ie -1) +0);
2161 nn = sdeg_mat(im) +1
2162
2163 allocate(ct(nn),ww(nn),dd(nn,nn))
2165 rho = prop_mat(im,1)
2166 lambda = prop_mat(im,2)
2167 mu = prop_mat(im,3)
2168
2169 if (n_test.gt.0) then
2170
2171 do ip = 1, n_test
2172
2173 fn = 0
2174 do ifun = 1,nfunc
2175 if (fun_test(ip).eq.tag_func(ifun)) fn = ifun
2176 enddo
2177
2178 if (fn.gt.0) then
2179
2180 rho = prop_mat(im,1); lambda = prop_mat(im,2); mu = prop_mat
2181
2182
2183 do k = 1,nn
2184 do j = 1,nn
2185 do i = 1,nn
2186 dxdx = alfa11(ie) +beta12(ie)*ct(k) &
2187 + beta13(ie)*ct(j) +gamma1(ie)*ct(j)*ct(k)
2188 dydx = alfa21(ie) +beta22(ie)*ct(k) &
2189 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
2190 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
2191 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
2192
2193 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
2194 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
2195 dydy = alfa22(ie) +beta21(ie)*ct(k) &
2196 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
2197 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
2198 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
2199
2200 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
2201 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
2202 dydz = alfa23(ie) +beta21(ie)*ct(j) &
2203 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
2204 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
2205 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
2206
2207 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
2208 - dydz * (dxdx*dzdy - dzdx*dxdy) &
2209 + dzdz * (dxdx*dydy - dydx*dxdy)
2210
2211 is = nn*nn*(k -1) +nn*(j -1) +i
2212 in = cs_loc(cs_loc(ie -1) +is)
2213
2214
2215 x = xs_loc(in); y = ys_loc(in); z = zs_loc(in
2216
2217
2218 fmat(fn,(3*(in -1) +1)) = fmat(fn,(3*(in -1)
2219 4.d0*pi**2.d0*dcos(pi*y)*dcos(pi*z)*dsin
2220 (2.d0*lambda - 8.d0*mu + 9.d0*rho - 4.d0
2221 + 8.d0*mu*dcos(pi*x)**2.d0 - 9.d0*rho*dcos
2222
2223 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1)
2224 4.d0*pi**2.d0*dcos(pi*x)*dcos(pi*z)*dsin
2225 (2.d0*lambda + 12.d0*mu - 9.d0*rho - 4.d0
2226 - 16.d0*mu*dcos(pi*y)**2.d0 + 9.d0*rho*dcos
2227
2228 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in -1)
2229 4.d0*pi**2.d0*dcos(pi*x)*dcos(pi*y)*dsin
2230 (2.d0*lambda + 12.d0*mu - 9.d0*rho - 4.d0
2231 - 16.d0*mu*dcos(pi*z)**2.d0 + 9.d0*rho*dcos
2232
2233
2234 enddo
2235 enddo
2236 enddo
2237 endif
2238 enddo
2239 endif
2240 deallocate(ct,ww,dd)
2241 enddo
2242
2243
2244if (cs_nnz_bc_loc.gt.0) then
2245nface_loc = cs_bc_loc(0) -1
2246
2247do ie = 1,nface_loc
2248
2249
2250nn = int(sqrt(real(cs_bc_loc(ie) - cs_bc_loc(ie -1) -1)))
2251
2252allocate(ct(nn),ww(nn),dd(nn,nn))
2253
2254
2255
2256ic1 = cs_bc_loc(cs_bc_loc(ie -1) +nn*nn)
2257ic2 = cs_bc_loc(cs_bc_loc(ie -1) +1)
2258ic3 = cs_bc_loc(cs_bc_loc(ie -1) +nn*(nn -1) +1)
2259ic4 = cs_bc_loc(cs_bc_loc(ie -1) +nn)
2260
2261l1x = 0.d0
2262l2x = 0.d0
2263l1y = 0.d0
2264l2y = 0.d0
2265l1z = 0.d0
2266l2z = 0.d0
2267area = 0.d0
2268
2269if (ic1 .ne. 0 .and. ic2 .ne. 0 .and. ic3 .ne. 0 .and. ic4 .ne. 0) then
2270
2271l1x = xs_loc(ic1) - xs_loc(ic2)
2272l1y = ys_loc(ic1) - ys_loc(ic2)
2273l1z = zs_loc(ic1) - zs_loc(ic2)
2274
2275l2x = xs_loc(ic3) - xs_loc(ic4)
2276l2y = ys_loc(ic3) - ys_loc(ic4)
2277l2z = zs_loc(ic3) - zs_loc(ic4)
2278
2279area = l1y*l2z -l1z*l2y +l1z*l2x -l1x*l2z +l1x*l2y -l1y*l2x
2280area = dabs(0.5d0*area)
2281
2282endif
2283
2284
2285if (nl_neux.gt.0) then
2286do il = 1,nl_neux
2287 if (tag_neux(il).eq.cs_bc_loc(cs_bc_loc(ie -1) +0)) then
2288 fn = 0
2289 do ifun = 1,nfunc
2290 if (fun_neux(il).eq.tag_func(ifun)) fn = ifun
2291 enddo
2292
2293 if (fn.gt.0) then
2294
2295 v1 = val_neux(il,1)
2296 v2 = val_neux(il,2)
2297 v3 = val_neux(il,3)
2298 v4 = val_neux(il,4)
2299
2300 do j = 1,nn
2301 do i = 1,nn
2302 is = nn*(j -1) +i
2303 in = cs_bc_loc(cs_bc_loc(ie -1) +is)
2304 v = 0.25*(1.0 -ct(i))*(1.0 -ct(j))*v1 &
2305 + 0.25*(1.0 +ct(i))*(1.0 -ct(j))*v2 &
2306 + 0.25*(1.0 +ct(i))*(1.0 +ct(j))*v3 &
2307 + 0.25*(1.0 -ct(i))*(1.0 +ct(j))*v4
2308
2309 term = 0.25*area*v * ww(i)*ww(j)
2310
2311 fmat(fn,(3*(in -1) +1)) = fmat(fn,(3*(in -1) +1)) + term
2312 enddo
2313 enddo
2314 endif
2315 endif
2316enddo
2317endif
2318
2319
2320if (nl_neuy.gt.0) then
2321do il = 1,nl_neuy
2322 if (tag_neuy(il).eq.cs_bc_loc(cs_bc_loc(ie -1) +0)) then
2323 fn = 0
2324 do ifun = 1,nfunc
2325 if (fun_neuy(il).eq.tag_func(ifun)) fn = ifun
2326 enddo
2327
2328 if (fn.gt.0) then
2329 v1 = val_neuy(il,1)
2330 v2 = val_neuy(il,2)
2331 v3 = val_neuy(il,3)
2332 v4 = val_neuy(il,4)
2333
2334 do j = 1,nn
2335 do i = 1,nn
2336 is = nn*(j -1) +i
2337 in = cs_bc_loc(cs_bc_loc(ie -1) +is)
2338 v = 0.25*(1.0 -ct(i))*(1.0 -ct(j))*v1 &
2339 + 0.25*(1.0 +ct(i))*(1.0 -ct(j))*v2 &
2340 + 0.25*(1.0 +ct(i))*(1.0 +ct(j))*v3 &
2341 + 0.25*(1.0 -ct(i))*(1.0 +ct(j))*v4
2342
2343 term = 0.25*area*v * ww(i)*ww(j)
2344
2345 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1) +2)) + term
2346 enddo
2347 enddo
2348 endif
2349 endif
2350enddo
2351endif
2352
2353
2354if (nl_neuz.gt.0) then
2355do il = 1,nl_neuz
2356 if (tag_neuz(il).eq.cs_bc_loc(cs_bc_loc(ie -1) +0)) then
2357 fn = 0
2358 do ifun = 1,nfunc
2359 if (fun_neuz(il).eq.tag_func(ifun)) fn = ifun
2360 enddo
2361
2362 if (fn.gt.0) then
2363 v1 = val_neuz(il,1)
2364 v2 = val_neuz(il,2)
2365 v3 = val_neuz(il,3)
2366 v4 = val_neuz(il,4)
2367
2368 do j = 1,nn
2369 do i = 1,nn
2370 is = nn*(j -1) +i
2371 in = cs_bc_loc(cs_bc_loc(ie -1) +is)
2372 v = 0.25*(1.0 -ct(i))*(1.0 -ct(j))*v1 &
2373 + 0.25*(1.0 +ct(i))*(1.0 -ct(j))*v2 &
2374 + 0.25*(1.0 +ct(i))*(1.0 +ct(j))*v3 &
2375 + 0.25*(1.0 -ct(i))*(1.0 +ct(j))*v4
2376
2377 term = 0.25*area*v * ww(i)*ww(j)
2378
2379 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in -1) +3)) + term
2380 enddo
2381 enddo
2382 endif
2383 endif
2384enddo
2385endif
2386
2387
2388
2389if (nl_neun.gt.0) then
2390do il = 1,nl_neun
2391 if (tag_neun(il).eq.cs_bc_loc(cs_bc_loc(ie -1) +0)) then
2392 fn = 0
2393 do ifun = 1,nfunc
2394 if (fun_neun(il).eq.tag_func(ifun)) fn = ifun
2395 enddo
2396
2397 if (fn.gt.0) then
2398
2399 rho = prop_mat(1,1); lambda = prop_mat(1,2); mu = prop_mat(1,3)
2400
2401 k1 = val_neun(il,1)
2402 k2 = val_neun(il,2)
2403 k3 = val_neun(il,3)
2404
2405 omega = val_neun(il,4)*2*pi
2406
2407 vp = dsqrt((lambda+2*mu)/rho)
2408 k3 = k3*omega/vp;
2409
2410
2411 v1 = 0;
2412 v2 = 0
2413 v3 = 0
2414
2415
2416 do j = 1,nn
2417 do i = 1,nn
2418 is = nn*(j -1) +i
2419 in = cs_bc_loc(cs_bc_loc(ie -1) +is)
2420
2421 x = xs_loc(in); y = ys_loc(in); z = zs_loc(in);
2422
2423 do index = 1,nelem_neun
2424 if (i4normal(index) .eq. ie) index_vector = index
2425 enddo
2426
2427 n1 = normal_nx_el_neun(index_vector)
2428 n2 = normal_ny_el_neun(index_vector)
2429 n3 = normal_nz_el_neun(index_vector)
2430
2431
2432
2433
2434 if (tag_func(fn) .eq. 202) then
2435 v1 = k3*lambda*n1*cos(k1*x + k2*y + k3*z) + k1*mu*n3*cos
2436 v2 = k3*lambda*n2*cos(k1*x + k2*y + k3*z) + k2*mu*n3*cos
2437 v3 = n3*(k3*lambda*cos(k1*x + k2*y + k3*z) + 2*k3*mu*cos
2438 + k1*mu*n1*cos(k1*x + k2*y + k3*z) + k2*mu*n2*cos(k1*x
2439
2440 v1 = -v1; v2 = -v2; v3 = -v3;
2441
2442 elseif(tag_func(fn) .eq. 203) then
2443 v1 = - k3*lambda*n1*sin(k1*x + k2*y + k3*z) - k1*mu*n3*sin
2444 v2 = - k3*lambda*n2*sin(k1*x + k2*y + k3*z) - k2*mu*n3*sin
2445 v3 = - n3*(k3*lambda*sin(k1*x + k2*y + k3*z) + 2*k3*mu*sin
2446 - k1*mu*n1*sin(k1*x + k2*y + k3*z) - k2*mu*n2*sin(k1
2447
2448 v1 = -v1; v2 = -v2; v3 = -v3;
2449 else
2450 write(*,*) 'Case not implemented!!! Check inputs.'
2451
2452 endif
2453
2454 if (n_test.gt.0) then
2455
2456 v1 = n1*(lambda*(2*pi*cos(pi*y)*sin(2*pi*x)*sin(pi*y)*sin(
2457 - 2*pi*cos(pi*x)*sin(pi*x)*sin(2*pi*y)*sin(2*pi*z) &
2458 + 2*pi*cos(pi*z)*sin(2*pi*x)*sin(2*pi*y)*sin(pi*z)) -
2459 4*mu*pi*cos(pi*x)*sin(pi*x)*sin(2*pi*y)*sin(2*pi*z))
2460 + mu*n2*(2*pi*cos(2*pi*x)*sin(pi*y)**2*sin(2*pi*z) &
2461 - 2*pi*cos(2*pi*y)*sin(pi*x)**2*sin(2*pi*z)) &
2462 + mu*n3*(2*pi*cos(2*pi*x)*sin(2*pi*y)*sin(pi*z)**2 &
2463 - 2*pi*cos(2*pi*z)*sin(pi*x)**2*sin(2*pi*y))
2464 v2 = n2*(lambda*(2*pi*cos(pi*y)*sin(2*pi*x)*sin(pi*y)*sin(2
2465 - 2*pi*cos(pi*x)*sin(pi*x)*sin(2*pi*y)*sin(2*pi*z) &
2466 + 2*pi*cos(pi*z)*sin(2*pi*x)*sin(2*pi*y)*sin(pi*z)) &
2467 + 4*mu*pi*cos(pi*y)*sin(2*pi*x)*sin(pi*y)*sin(2*pi*z)) &
2468 + mu*n1*(2*pi*cos(2*pi*x)*sin(pi*y)**2*sin(2*pi*z) &
2469 - 2*pi*cos(2*pi*y)*sin(pi*x)**2*sin(2*pi*z)) &
2470 + mu*n3*(2*pi*cos(2*pi*y)*sin(2*pi*x)*sin(pi*z)**2 &
2471 + 2*pi*cos(2*pi*z)*sin(2*pi*x)*sin(pi*y)**2)
2472 v3 = n3*(lambda*(2*pi*cos(pi*y)*sin(2*pi*x)*sin(pi*y)*sin(2
2473 - 2*pi*cos(pi*x)*sin(pi*x)*sin(2*pi*y)*sin(2*pi*z) &
2474 + 2*pi*cos(pi*z)*sin(2*pi*x)*sin(2*pi*y)*sin(pi*z)) &
2475 + 4*mu*pi*cos(pi*z)*sin(2*pi*x)*sin(2*pi*y)*sin(pi*z)) &
2476 + mu*n1*(2*pi*cos(2*pi*x)*sin(2*pi*y)*sin(pi*z)**2 &
2477 - 2*pi*cos(2*pi*z)*sin(pi*x)**2*sin(2*pi*y)) &
2478 + mu*n2*(2*pi*cos(2*pi*y)*sin(2*pi*x)*sin(pi*z)**2 &
2479 + 2*pi*cos(2*pi*z)*sin(2*pi*x)*sin(pi*y)**2)
2480 endif
2481
2482 fmat(fn,(3*(in -1) +1)) = fmat(fn,(3*(in -1) +1)) + 0.2
2483 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1) +2)) + 0.2
2484 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in -1) +3)) + 0.2
2485
2486 enddo
2487 enddo
2488 endif
2489 endif
2490enddo
2491
2492endif
2493
2494
2495
2496if (nl_dirx.gt.0) then
2497do il = 1,nl_dirx
2498 if (tag_dirx(il).eq.cs_bc_loc(cs_bc_loc(ie -1) +0)) then
2499 fn = 0
2500 do ifun = 1,nfunc
2501 if (fun_dirx(il).eq.tag_func(ifun)) fn = ifun
2502 enddo
2503
2504 if (fn.gt.0) then
2505 v1 = val_dirx(il,1)
2506 v2 = val_dirx(il,2)
2507 v3 = val_dirx(il,3)
2508 v4 = val_dirx(il,4)
2509
2510 do j = 1,nn
2511 do i = 1,nn
2512 is = nn*(j -1) +i
2513 in = cs_bc_loc(cs_bc_loc(ie -1) +is)
2514 v = 0.25*(1.0 -ct(i))*(1.0 -ct(j))*v1 &
2515 + 0.25*(1.0 +ct(i))*(1.0 -ct(j))*v2 &
2516 + 0.25*(1.0 +ct(i))*(1.0 +ct(j))*v3 &
2517 + 0.25*(1.0 -ct(i))*(1.0 +ct(j))*v4
2518
2519 fmat(fn,(3*(in -1) +1)) = v
2520 enddo
2521 enddo
2522 endif
2523 endif
2524enddo
2525endif
2526
2527if (nl_diry.gt.0) then
2528do il = 1,nl_diry
2529 if (tag_diry(il).eq.cs_bc_loc(cs_bc_loc(ie -1) +0)) then
2530 fn = 0
2531 do ifun = 1,nfunc
2532 if (fun_diry(il).eq.tag_func(ifun)) fn = ifun
2533 enddo
2534
2535 if (fn.gt.0) then
2536 v1 = val_diry(il,1)
2537 v2 = val_diry(il,2)
2538 v3 = val_diry(il,3)
2539 v4 = val_diry(il,4)
2540
2541 do j = 1,nn
2542 do i = 1,nn
2543 is = nn*(j -1) +i
2544 in = cs_bc_loc(cs_bc_loc(ie -1) +is)
2545 v = 0.25*(1.0 -ct(i))*(1.0 -ct(j))*v1 &
2546 + 0.25*(1.0 +ct(i))*(1.0 -ct(j))*v2 &
2547 + 0.25*(1.0 +ct(i))*(1.0 +ct(j))*v3 &
2548 + 0.25*(1.0 -ct(i))*(1.0 +ct(j))*v4
2549
2550 fmat(fn,(3*(in -1) +2)) = v
2551 enddo
2552 enddo
2553 endif
2554 endif
2555enddo
2556endif
2557
2558
2559if (nl_dirz.gt.0) then
2560do il = 1,nl_dirz
2561 if (tag_dirz(il).eq.cs_bc_loc(cs_bc_loc(ie -1) +0)) then
2562 fn = 0
2563 do ifun = 1,nfunc
2564 if (fun_dirz(il).eq.tag_func(ifun)) fn = ifun
2565 enddo
2566
2567 if (fn.gt.0) then
2568 v1 = val_dirz(il,1)
2569 v2 = val_dirz(il,2)
2570 v3 = val_dirz(il,3)
2571 v4 = val_dirz(il,4)
2572
2573 do j = 1,nn
2574 do i = 1,nn
2575 is = nn*(j -1) +i
2576 in = cs_bc_loc(cs_bc_loc(ie -1) +is)
2577 v = 0.25*(1.0 -ct(i))*(1.0 -ct(j))*v1 &
2578 + 0.25*(1.0 +ct(i))*(1.0 -ct(j))*v2 &
2579 + 0.25*(1.0 +ct(i))*(1.0 +ct(j))*v3 &
2580 + 0.25*(1.0 -ct(i))*(1.0 +ct(j))*v4
2581
2582 fmat(fn,(3*(in -1) +3)) = v
2583 enddo
2584 enddo
2585 endif
2586 endif
2587enddo
2588endif
2589deallocate(ct,dd,ww)
2590enddo
2591endif
2592
2593if (nl_neun.gt.0) deallocate(i4normal, normal_nx_el_neun, normal_ny_el_neun
2594if (nl_poiz.gt.0) deallocate(node_counter_poiz)
2595if (nl_poiy.gt.0) deallocate(node_counter_poiy)
2596if (nl_poix.gt.0) deallocate(node_counter_poix)
2597
2598
2599
2600
2601
2602
2603
2604
2605
2606
2607
2608
2609return
subroutine get_elem_from_surf(nnz_bc, cs_bc, n1, n2, n3, n4, ie_surf)
Find element index from verteces number.
subroutine get_indloc_from_indglo(local_el, nel_loc, ie, ic)
Returns local id from global id.
subroutine get_mech_prop_nh_enhanced(ie, nn, nn_loc, cs_nnz_loc, cs_loc, rho_nhe, lambda_nhe, mu_nhe, qs_nh, fmax, rho_el, lambda_el, mu_el, gamma_el)
...Not-Honoring Enhanced (NHE) Implementation
subroutine get_minvalues(i_glo, v_glo, n_glo, i_loc, n_loc, np)
Computes positions of minimum values.
subroutine get_nearest_node(n, xs, ys, zs, xt, yt, zt, nt, dist_min)
Computes the nearest node with respet to (xt,yt,zt).
subroutine get_node_from_face(nb_nodes, nb_nz, con_spc, nb_load, l
Computes total number of boundary nodes and second derivative on a given point x.
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
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.
real *8 function phi(ng, csii, csij)
Evaluates Lengendre basis function on a specific point.
integer, parameter exit_surf_notfound