SPEED
MAKE_EXTINT_FORCES.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
93
94
95 subroutine make_extint_forces(nnod_loc,xs_loc,ys_loc,zs_loc,local_n_num,cs_nnz_loc,cs_loc,&
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, & !val_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
142 use speed_exit_codes
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
154 integer*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
156 integer*4 :: nl_sism, srcmodflag, szsism
157 integer*4 :: nl_expl
158 integer*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
164 integer*4 :: nhexa
165 integer*4 :: sum_node_bottom,sum_node_bottom1,sum_node_bottom2
166 integer*4 :: sum_node_bottom_first3
167 integer*4 :: last_node_bottom,last_node_bottom1,last_node_bottom2
168 integer*4 :: C,sit
169 integer*4 :: isism,ipsism
170 integer*4 :: length_cns, sum_length_cns
171 integer*4 :: max_num_ns
172 integer*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
178 integer*4 :: nnode_neuN,nelem_neuN
179 integer*4 :: ne1,ne2,ne3,ne4
180 integer*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
188 integer*4, dimension(nl_expl) :: num_ne
189 integer*4, dimension(6) :: num_sit
190 integer*4, dimension(nl_sism) :: num_ns
191 integer*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
213 integer*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
235 integer*4, dimension(max_num_ns,nl_sism) :: sour_ns
236 integer*4, dimension(max_num_ne,nl_expl) :: sour_ne
237
238 integer*4, dimension(:), allocatable :: node_counter_poiX,node_counter_poiY,&
239 node_counter_poiz, recv_poiz, &
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
247 real*8 :: dist, phi
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, k2, k3
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,z_tra_glo
270 real*8, dimension(:), allocatable :: xt_tra,yt_tra,zt_tra,xt_tra_glo,yt_tra_glo,zt_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, gamma_el
275 real*8, dimension(nm) :: tref_mat
276 real*8, dimension(ne_loc) :: alfa11,alfa12,alfa13, alfa21,alfa22,alfa23, alfa31,alfa32,alfa33
277 real*8, dimension(ne_loc) :: beta11,beta12,beta13, beta21,beta22,beta23, beta31,beta32,beta33
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 !real*8, dimension(n_test,10) :: val_fun_test
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! Initialization "seismic moment" variables - begin
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! Initialazation "seismic moment" variables - end
335
336! Initialization "explosive source" variables - begin
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! Initialazation "explosive source" variables - end
349
350
351
352!*************************************************************************************
353! Calculation of the normal vector with respect to the surface element in which the
354! Neumann load is enforced
355! normal_nx_el_neuN, normal_ny_el_neuN, normal_nz_el_neuN = are the 3 components
356! of the normal vector
357!*************************************************************************************
358
359 if (nl_neun.gt.0) then
360
361
362 nelem_neun = 0
363 allocate(i4count(nnod_loc))
364 i4count = 0
365
366 call get_node_from_face(nnod_loc,cs_nnz_bc_loc,cs_bc_loc,nl_neun,tag_neun,&
367 nnode_neun, i4count, local_n_num)
368
369
370
371 ! nnode_neuN = node number in which the Neumann load is applied
372 ! i4count(i) different from 0 if on the i-th the Neumann load is applied
373 ! ne_loc = number of local elements
374 ! ielem = global index for the element
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 -1) +1)
383 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
384 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
385 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
386
387 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
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 -1) +1)
392 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
393 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
394 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
395
396 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
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 -1) +1)
401 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
402 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
403 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
404
405 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
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 -1) +1)
410 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
411 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
412 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
413
414 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
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 -1) +1)
419 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
420 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
421 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
422
423 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
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 -1) +1)
428 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
429 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
430 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
431
432 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
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) +1)
459 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
460 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
461 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
462
463 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
464
465 !face x = -1
466
467 nelem_neun = nelem_neun +1
468
469 call get_elem_from_surf(cs_nnz_bc_loc, cs_bc_loc, ne1, ne2, ne3, ne4, ie_surf)
470 if(ie_surf .eq. 0) then
471 write(*,*) 'surf not found'
472 call exit(exit_surf_notfound)
473 endif
474 i4normal(nelem_neun) = ie_surf
475
476 call make_normal(1,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3), xs_loc(ne4), &
477 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3), ys_loc(ne4), &
478 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3), zs_loc(ne4) ,&
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 -1) +1)
490 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
491 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
492 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
493
494 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
495
496 ! face y = -1
497 nelem_neun = nelem_neun +1
498
499 call get_elem_from_surf(cs_nnz_bc_loc,cs_bc_loc, ne1, ne2, ne3, ne4, ie_surf)
500 if(ie_surf .eq. 0) then
501 write(*,*) 'surf not found'
502 call exit(exit_surf_notfound)
503 endif
504 i4normal(nelem_neun) = ie_surf
505
506 call make_normal(2,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3), xs_loc(ne4), &
507 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3), ys_loc(ne4), &
508 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3), zs_loc(ne4) ,&
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 -1) +1)
519 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
520 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
521 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
522
523 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
524
525 !face z = -1
526 nelem_neun = nelem_neun +1
527
528 call get_elem_from_surf(cs_nnz_bc_loc,cs_bc_loc, ne1, ne2, ne3, ne4, ie_surf)
529 if(ie_surf .eq. 0) then
530 write(*,*) 'surf not found'
531 call exit(exit_surf_notfound)
532 endif
533 i4normal(nelem_neun) = ie_surf
534
535 call make_normal(3,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3), xs_loc(ne4), &
536 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3), ys_loc(ne4), &
537 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3), zs_loc(ne4) ,&
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 -1) +1)
548 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
549 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
550 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
551
552 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
553
554 ! face x = 1
555
556 nelem_neun = nelem_neun +1
557
558 call get_elem_from_surf(cs_nnz_bc_loc,cs_bc_loc, ne1, ne2, ne3, ne4, ie_surf)
559 if(ie_surf .eq. 0) then
560 write(*,*) 'surf not found'
561 call exit(exit_surf_notfound)
562 endif
563
564 i4normal(nelem_neun) = ie_surf
565
566 call make_normal(4,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3), xs_loc(ne4), &
567 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3), ys_loc(ne4), &
568 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3), zs_loc(ne4) ,&
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 -1) +1)
579 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
580 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
581 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
582
583 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
584
585 ! face y = 1
586
587 nelem_neun = nelem_neun +1
588
589 call get_elem_from_surf(cs_nnz_bc_loc,cs_bc_loc, ne1, ne2, ne3, ne4, ie_surf)
590 if(ie_surf .eq. 0) then
591 write(*,*) 'surf not found'
592 call exit(exit_surf_notfound)
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), xs_loc(ne4), &
599 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3), ys_loc(ne4), &
600 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3), zs_loc(ne4) ,&
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 -1) +1)
611 ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
612 ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
613 ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
614
615 if ((i4count(ne1).ne.0).and.(i4count(ne2).ne.0) .and.(i4count(ne3).ne.0) .and.(i4count(ne4).ne.0)) then
616
617 ! face z = 1
618 nelem_neun = nelem_neun +1
619
620 call get_elem_from_surf(cs_nnz_bc_loc,cs_bc_loc, ne1, ne2, ne3, ne4, ie_surf)
621 if(ie_surf .eq. 0) then
622 write(*,*) 'surf not found'
623 call exit(exit_surf_notfound)
624 endif
625 i4normal(nelem_neun) = ie_surf
626
627 call make_normal(6,xs_loc(ne1), xs_loc(ne2), xs_loc(ne3), xs_loc(ne4), &
628 ys_loc(ne1), ys_loc(ne2), ys_loc(ne3), ys_loc(ne4), &
629 zs_loc(ne1), zs_loc(ne2), zs_loc(ne3), zs_loc(ne4) ,&
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 !if(nl_neuN.gt.0)
646
647 ! *******************************************************************************************
648
649 ! Fmat(nfunc, 3*nnode_dom) ---> Fmat(nfunc, nnod_loc)
650
651 fmat = 0.0d0
652
653 !*********************************************************************************************
654 ! If nl_poiX is different from 0 first compute the nearest local nodes to the one in which the
655 ! load is applied and then compute the nearest global node by comparing the resutls on each
656 ! mpi process
657
658 if (nl_poix .gt. 0) then
659
660 allocate(val_locx(nl_poix), val_glox(nl_poix*mpi_np), ind_locx(nl_poix), ind_glox(nl_poix*mpi_np))
661
662 do i = 1,nl_poix
663 call get_nearest_node(nnod_loc,xs_loc,ys_loc,zs_loc, val_poix(i,1),val_poix(i,2),val_poix(i,3),node_poix(i), dist)
664 val_locx(i) = dist
665 ind_locx(i) = local_n_num(node_poix(i)) !local_n_num(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, nl_poix, speed_double, mpi_comm, mpi_ierr)
672 call mpi_allgather(ind_locx, nl_poix, speed_integer, ind_glox, nl_poix, speed_integer, mpi_comm, mpi_ierr)
673
674 ind_locx = 0
675
676 call get_minvalues(ind_glox, val_glox, nl_poix*mpi_np, ind_locx, nl_poix, mpi_np)
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) + 1
701 endif
702 enddo
703 enddo
704 enddo
705 enddo
706
707 enddo
708
709
710 ! write(*,*) mpi_id, node_counter_poiX
711 call mpi_barrier(mpi_comm,mpi_ierr)
712 call mpi_allreduce(node_counter_poix,recv_poix, nl_poix, speed_integer, mpi_sum, mpi_comm, mpi_ierr)
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(nl_poiy), ind_gloy(nl_poiy*mpi_np))
722
723 do i = 1,nl_poiy
724 call get_nearest_node(nnod_loc,xs_loc,ys_loc,zs_loc, val_poiy(i,1),val_poiy(i,2),val_poiy(i,3),node_poiy(i), dist)
725 val_locy(i) = dist
726 ind_locy(i) = local_n_num(node_poiy(i)) !local_n_num(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, nl_poiy, speed_double, mpi_comm, mpi_ierr)
732 call mpi_allgather(ind_locy, nl_poiy, speed_integer, ind_gloy, nl_poiy, speed_integer, mpi_comm, mpi_ierr)
733
734 call get_minvalues(ind_gloy,val_gloy, nl_poiy*mpi_np, ind_locy, nl_poiy, mpi_np)
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) + 1
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, mpi_sum, mpi_comm, mpi_ierr)
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(nl_poiz), ind_gloz(nl_poiz*mpi_np))
781
782 do i = 1,nl_poiz
783 call get_nearest_node(nnod_loc,xs_loc,ys_loc,zs_loc, val_poiz(i,1),val_poiz(i,2),val_poiz(i,3),node_poiz(i), dist)
784 val_locz(i) = dist
785 ind_locz(i) = local_n_num(node_poiz(i)) !local_n_num(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, nl_poiz, speed_double, mpi_comm, mpi_ierr)
791 call mpi_allgather(ind_locz, nl_poiz, speed_integer, ind_gloz, nl_poiz, speed_integer, mpi_comm, mpi_ierr)
792 call get_minvalues(ind_gloz, val_gloz, nl_poiz*mpi_np, ind_locz, nl_poiz, mpi_np)
793
794 do i = 1, nl_poiz
795 node_poiz(i) = ind_gloz(ind_locz(i))
796 ! write(*,*) i, ind_locZ(i), ind_gloZ(ind_locZ(i)), node_poiZ(i)
797 enddo
798 ! write(*,*) node_poiZ
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) + 1
819 endif
820 enddo
821 enddo
822 enddo
823 enddo
824
825 enddo
826
827
828 !write(*,*) mpi_id, node_counter_poiZ
829 call mpi_barrier(mpi_comm,mpi_ierr)
830 call mpi_allreduce(node_counter_poiz, recv_poiz, nl_poiz, speed_integer, mpi_sum, mpi_comm, mpi_ierr)
831 node_counter_poiz = recv_poiz;
832 deallocate(recv_poiz);
833 !write(*,*) mpi_id, node_counter_poiZ
834 call mpi_barrier(mpi_comm,mpi_ierr)
835
836 !stop
837 endif
838
839 ne_loc = cs_loc(0) -1
840
841
842
843 do ie = 1, ne_loc
844 !write(*,*) ie, '/', ne_loc
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))
849 call make_lgl_nw(nn,ct,ww,dd)
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,nn,nn), gamma_el(nn,nn,nn))
857 call get_mech_prop_nh_enhanced(ie, nn, nnod_loc, cs_nnz_loc, cs_loc, &
858 rho_nhe, lambda_nhe, mu_nhe, 100.d0, 3.d0, &
859 rho_el, lambda_el, mu_el, gamma_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 ! Point load X
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
875 call get_indloc_from_indglo(local_n_num, nnod_loc, node_poix(ip), ic1)
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(ip)) then
886 fmat(fn,(3*(in -1) +1)) = fmat(fn,(3*(in -1) +1)) + val_poix(ip,4)/node_counter_poix(ip)
887 ! node_counter_poiX = node_counter_poiX + 1;
888
889 ! SDOF SS (Need to define, building ids corresponding to each node)
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 ! Point load Y
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
909 call get_indloc_from_indglo(local_n_num, nnod_loc, node_poiy(ip), ic1)
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(ip)) then
919 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1) +2)) + val_poiy(ip,4)/node_counter_poiy(ip)
920 endif
921 enddo
922 enddo
923 enddo
924 endif
925 endif
926 enddo
927 endif !(nl_poiY .gt. 0)
928
929 if (nl_poiz.gt.0) then ! Point load Z
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
938 call get_indloc_from_indglo(local_n_num, nnod_loc, node_poiz(ip), ic1)
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 -1) +3)) &
951 + val_poiz(ip,4)/node_counter_poiz(ip)
952
953 endif
954
955 enddo
956 enddo
957 enddo
958 endif
959 endif
960 enddo
961 endif !(nl_poiZ .gt. 0)
962
963
964 !#############################################################################
965 !################## PLANE WAVE BEGIN ##############
966 !#############################################################################
967
968 if (nl_plax.gt.0) then ! Plane Wave X
969
970 do ipl = 1,nl_plax
971
972 if (tag_mat(im).eq.tag_plax(ipl)) then !Check on the Plane Wave Material - start
973
974 ! Recognizing bottom face - begin
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
983 call get_indloc_from_indglo(local_n_num, nnod_loc, con_hexa(ie_glob,2), ic1)
984 call get_indloc_from_indglo(local_n_num, nnod_loc, con_hexa(ie_glob,i+1), ic2)
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 ! Recognizing bottom face - end
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 !Table of the possible situation:
1021 !hypothesis:
1022 !- the plane wave load source is placed on a rectangular element, aligned
1023 ! along the principal axis (x,y,z). It does NOT work properly on element
1024 ! with different height (z co-ordinate MUST BE EQUAL)!!!
1025 !- the node orientation is ALWAYS counter-clockwise!!!
1026 !
1027 ! Counter Clock wise means: | Looking from the top:
1028 ! | The number of nodes increses
1029 ! ^ Z | going from the x axis versus y axis
1030 ! | | 1 - 2 - 3 - 4
1031 ! | |
1032 ! 5------8 |
1033 ! /| /| | Z
1034 ! / | / | | o-----------------> Y
1035 ! 6------7 | | |
1036 ! | 1---|--4 ------> Y | | 1-------4
1037 ! | / | / | | | |
1038 ! |/ |/ | | | |
1039 ! 2------3 | | | |
1040 ! / | | 2-------3
1041 ! / | |
1042 ! / | v X
1043 ! v X |
1044 !
1045 !
1046 ! Possible position of the BOTTOM FACE
1047 !
1048 ! S1 S2 S3
1049 !
1050 ! 5------8 8------4 4------2
1051 ! /| /| /| /| /------/|
1052 ! / | / | / | //| /------/ |
1053 ! 6------7 | 7------3//| 3------1 |
1054 ! | 1---|--4 | 5---|//2 | 8---|--5
1055 ! | /----|-/ | / |// | / | /
1056 ! |/-----|/ |/ |/ |/ |/
1057 ! 2------3 6------1 7------6
1058 !
1059 ! sum_node_bottom
1060 ! 1+2+3+4=10 5+6+1+2=14 8+7+6+5=26
1061 !
1062 ! S4 S5 S6
1063 !
1064 ! 2------5 6------7 1------4
1065 ! /| /| /| /| /|-----/|
1066 ! //| / | / | / | / |----/-|
1067 ! 1------6 | 2------3 | 5------8--|
1068 ! |//4---|--8 |- 5---|--8 | 2---|--3
1069 ! |// | / |------| / | / | /
1070 ! |/ |/ |------|/ |/ |/
1071 ! 3------7 1------4 6------7
1072 !
1073 !
1074 ! 4+3+7+8=22 5+1+4+8=18 2+6+7+3=18
1075 !
1076
1077 !****** Situation S1 or S3 - begin ******
1078
1079 if ((sum_node_bottom.eq.10).or.(sum_node_bottom.eq.26)) then
1080 !
1081 !*** Situation S1 or S3 ***
1082 !
1083 ! e.g.: spectral degree = 3, nn = 4
1084 ! [#] = macro nodes
1085 ! # = micro nodes
1086 !
1087 ! i = 1,nn = 1,4
1088 ! j = 1,nn = 1,4
1089 ! z = ((nn-1)/2)+1 + num_sit(sit) = 2
1090 !
1091 !
1092 ! *******************************************************
1093 ! EXAMPLE DRAWING STILL ON WORKING!!!
1094 !
1095 !
1096 ! ^ Z
1097 ! |
1098 ! |
1099 ! 5------8
1100 ! /| /|
1101 ! / | / |
1102 ! 6------7 |
1103 ! |[1]
1104 ! 1---|--4 ------> Y
1105 ! |/ | /
1106 ! 2
1107 ! /
1108 ! 3
1109 ! /
1110 ! 4
1111 ! [2]-- 8--12--16[3]
1112 ! 2------3
1113 ! /
1114 ! /
1115 ! /
1116 ! v X
1117 !
1118 !
1119 ! S1
1120 ! [2]4--3--2--1[1]
1121 ! | |
1122 ! |8---7--6---5| Applied plane wave load on this surface
1123 ! | |
1124 ! |12--11-10--9|
1125 ! | |
1126 ! [3]16-15-14-13[4]
1127 !
1128 ! e.g.: spectral degree = 4, nn = 5
1129 !
1130 ! S1
1131 ! [2]5--4--3--2--1[1]
1132 ! | |
1133 ! |10--9--8--7---6|
1134 ! | |
1135 ! ->|15-14-13-12--11|<- Applied plane wave load
1136 ! | |
1137 ! |20-19-18-17--16|
1138 ! | |
1139 ! [3]25-24-23-22-21[4]
1140 !
1141 ! *******************************************************
1142 ! EXAMPLE DRAWING STILL ON WORKING!!!
1143 !
1144 !
1145
1146 !*** Situation S1 ***
1147
1148 if (sum_node_bottom.eq.10) then
1149
1150 sit = 1
1151
1152 !*** Situation S3 ***
1153
1154 else
1155
1156 sit = 2
1157
1158 endif
1159
1160 ! Even spectral polynomial degree
1161 if (mod(nn-1,2).eq.0) then
1162 k = ((nn-1)/2)+1
1163
1164 ! Odd spectral polynomial degree
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(ie)*ct(j)*ct(k)
1174 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1175 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1176 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1177 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1178 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1179 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1180 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1181 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1182 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1183 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1184 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1185 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1186 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1187 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1188 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1189 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
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(i) * ww(j) * dabs(det_j) / ellez
1197 fmat(fn,(3*(in -1) +1)) = fmat(fn,(3*(in -1) +1)) + term
1198
1199 prova = 3*(in -1) +1
1200
1201 enddo
1202 enddo
1203
1204
1205
1206
1207 !****** Situation S2 or S4 - begin ******
1208
1209 elseif ((sum_node_bottom.eq.14).or.(sum_node_bottom.eq.22)) then
1210
1211 !
1212 !*** Situation S2 ***
1213
1214 if (sum_node_bottom.eq.14) then
1215
1216 sit = 1
1217
1218 !*** Situation S4 ***
1219
1220 else
1221
1222 sit = 2
1223
1224 endif
1225
1226 ! Even spectral polynomial degree
1227 if (mod(nn-1,2).eq.0) then
1228 j = ((nn-1)/2)+1
1229
1230 ! Odd spectral polynomial degree
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)*ct(j)*ct(k)
1239 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1240 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1241 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1242 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1243 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1244 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1245 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1246 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1247 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1248 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1249 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1250 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1251 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1252 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1253 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1254 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
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(i) * ww(k) * dabs(det_j) / ellez
1263
1264 prova = 3*(in -1) +1
1265
1266
1267 fmat(fn,(3*(in -1) +1)) = fmat(fn,(3*(in -1) +1)) + term
1268 enddo
1269 enddo
1270
1271
1272
1273 !****** Situation S5 or S6 - begin ******
1274
1275 elseif (sum_node_bottom.eq.18) then
1276
1277 !*** Situation S5 ***
1278
1279 if ((sum_node_bottom_first3.eq.10) .or.(sum_node_bottom_first3.eq.13) &
1280 .or.(sum_node_bottom_first3.eq.14)) then
1281
1282 sit = 1
1283
1284 !*** Situation S6 ***
1285
1286 else
1287
1288 sit = 2
1289
1290 endif
1291
1292 ! Even spectral polynomial degree
1293 if (mod(nn-1,2).eq.0) then
1294 i = ((nn-1)/2)+1
1295
1296 ! Odd spectral polynomial degree
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)*ct(j)*ct(k)
1305 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1306 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1307 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1308 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1309 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1310 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1311 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1312 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1313 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1314 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1315 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1316 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1317 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1318 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1319 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1320 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
1321 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
1322 - dydz * (dxdx*dzdy - dzdx*dxdy) &
1323 + dzdz * (dxdx*dydy - dydx*dxdy)
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(j) * ww(k) * dabs(det_j) / ellez
1329 fmat(fn,(3*(in -1) +1)) = fmat(fn,(3*(in -1) +1)) + term
1330
1331
1332 enddo
1333 enddo
1334
1335 endif
1336
1337
1338 endif !if (fn.gt.0)
1339 endif !Check on the Plane Wave Material - end
1340
1341 enddo !ipl = 1,nl_plaX
1342 endif ! Plane Wave X - end
1343
1344 if (nl_play.gt.0) then ! Plane Wave Y
1345 do ipl = 1,nl_play
1346
1347 if (tag_mat(im).eq.tag_play(ipl)) then !Check on the Plane Wave Material - start
1348
1349 ! Recognizing bottom face - begin
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)
1357 call get_indloc_from_indglo(local_n_num, nnod_loc, con_hexa(ie_glob,2), ic1)
1358 call get_indloc_from_indglo(local_n_num, nnod_loc, con_hexa(ie_glob,i+1), ic2)
1359
1360 if(ic1 .ne. 0 .and. ic2 .ne. 0) then
1361 if (dabs(zs_loc(ic1) - zs_loc(ic2)) .le. 1.d-4) then
1362 !if (zs_loc(ic1) .eq. zs_loc(ic2)) then
1363 last_node_bottom1 = i
1364 sum_node_bottom1 = sum_node_bottom1 + i
1365 z_node_bottom1 = zs_loc(ic1)
1366 else
1367 last_node_bottom2 = i
1368 sum_node_bottom2 = sum_node_bottom2 + i
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 - last_node_bottom1
1377 else
1378 sum_node_bottom = sum_node_bottom2
1379 sum_node_bottom_first3 = sum_node_bottom2 - last_node_bottom2
1380 endif
1381
1382
1383 ! Recognizing bottom face - end
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 !****** Situation S1 or S3 - begin ******
1396
1397 if ((sum_node_bottom.eq.10).or.(sum_node_bottom.eq.26)) then
1398
1399
1400 !*** Situation S1 ***
1401
1402 if (sum_node_bottom.eq.10) then
1403
1404 sit = 1
1405
1406 !*** Situation S3 ***
1407
1408 else
1409
1410 sit = 2
1411
1412 endif
1413
1414 ! Even spectral polynomial degree
1415 if (mod(nn-1,2).eq.0) then
1416 k = ((nn-1)/2)+1
1417
1418 ! Odd spectral polynomial degree
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)*ct(j)*ct(k)
1427 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1428 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1429 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1430 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1431 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1432 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1433 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1434 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1435 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1436 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1437 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1438 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1439 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1440 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1441 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1442 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
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(i) * ww(j) * dabs(det_j) / ellez
1452 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1) +2)) + term
1453 enddo
1454 enddo
1455
1456
1457
1458
1459 !****** Situation S2 or S4 - begin ******
1460
1461 elseif ((sum_node_bottom.eq.14).or.(sum_node_bottom.eq.22)) then
1462
1463 !*** Situation S2 ***
1464
1465 if (sum_node_bottom.eq.14) then
1466
1467 sit = 1
1468
1469 !*** Situation S4 ***
1470
1471 else
1472
1473 sit = 2
1474
1475 endif
1476
1477 ! Even spectral polynomial degree
1478 if (mod(nn-1,2).eq.0) then
1479 j = ((nn-1)/2)+1
1480
1481 ! Odd spectral polynomial degree
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(j)*ct(k)
1491 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1492 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1493 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1494 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1495 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1496 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1497 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1498 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1499 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1500 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1501 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1502 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1503 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1504 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1505 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1506 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
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) * ww(k) * dabs(det_j) / ellez
1515 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1) +2)) + term
1516 enddo
1517 enddo
1518
1519
1520
1521 !****** Situation S5 or S6 - begin ******
1522
1523 elseif (sum_node_bottom.eq.18) then
1524
1525 !*** Situation S5 ***
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 !*** Situation S6 ***
1534
1535 else
1536
1537 sit = 2
1538
1539 endif
1540
1541 ! Even spectral polynomial degree
1542 if (mod(nn-1,2).eq.0) then
1543 i = ((nn-1)/2)+1
1544
1545 ! Odd spectral polynomial degree
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(j)*ct(k)
1555 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1556 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1557 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1558 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1559 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1560 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1561 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1562 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1563 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1564 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1565 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1566 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1567 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1568 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1569 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1570 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
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) * ww(k) * dabs(det_j) / ellez
1579 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1) +2)) + term
1580
1581 enddo
1582 enddo
1583
1584 endif
1585
1586
1587 endif
1588 endif !Check on the Plane Wave Material - end
1589
1590 enddo
1591 endif ! Plane Wave Y - end
1592
1593 if (nl_plaz.gt.0) then ! Plane Wave Z
1594 do ipl = 1,nl_plaz
1595
1596 if (tag_mat(im).eq.tag_plaz(ipl)) then !Check on the Plane Wave Material - start
1597
1598 ! Recognizing bottom face - begin
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)
1606 call get_indloc_from_indglo(local_n_num, nnod_loc, con_hexa(ie_glob,2), ic1)
1607 call get_indloc_from_indglo(local_n_num, nnod_loc, con_hexa(ie_glob,i+1), ic2)
1608
1609 if(ic1 .ne. 0 .and. ic2 .ne. 0) then
1610 if (dabs(zs_loc(ic1) - zs_loc(ic2)) .le. 1.d-4) then
1611 last_node_bottom1 = i
1612 sum_node_bottom1 = sum_node_bottom1 + i
1613 z_node_bottom1 = zs_loc(ic1)
1614 else
1615 last_node_bottom2 = i
1616 sum_node_bottom2 = sum_node_bottom2 + i
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 - last_node_bottom1
1625 else
1626 sum_node_bottom = sum_node_bottom2
1627 sum_node_bottom_first3 = sum_node_bottom2 - last_node_bottom2
1628 endif
1629
1630
1631 ! Recognizing bottom face - end
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 !****** Situation S1 or S3 - begin ******
1644
1645 if ((sum_node_bottom.eq.10).or.(sum_node_bottom.eq.26)) then
1646
1647
1648 !*** Situation S1 ***
1649
1650 if (sum_node_bottom.eq.10) then
1651
1652 sit = 1
1653
1654 !*** Situation S3 ***
1655
1656 else
1657
1658 sit = 2
1659
1660 endif
1661
1662 ! Even spectral polynomial degree
1663 if (mod(nn-1,2).eq.0) then
1664 k = ((nn-1)/2)+1
1665
1666 ! Odd spectral polynomial degree
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(k) &
1674 + beta13(ie)*ct(j) +gamma1(ie)*ct(j)*ct(k)
1675 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1676 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1677 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1678 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1679 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1680 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1681 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1682 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1683 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1684 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1685 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1686 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1687 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1688 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1689 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1690 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
1691 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
1692 - dydz * (dxdx*dzdy - dzdx*dxdy) &
1693 + dzdz * (dxdx*dydy - dydx*dxdy)
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) * ww(i) * ww(j) * dabs(det_j) / ellez
1699 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in -1) +3)) + term
1700 enddo
1701 enddo
1702
1703
1704
1705
1706 !****** Situation S2 or S4 - begin ******
1707
1708 elseif ((sum_node_bottom.eq.14).or.(sum_node_bottom.eq.22)) then
1709
1710 !*** Situation S2 ***
1711
1712 if (sum_node_bottom.eq.14) then
1713
1714 sit = 1
1715
1716 !*** Situation S4 ***
1717
1718 else
1719
1720 sit = 2
1721
1722 endif
1723
1724 ! Even spectral polynomial degree
1725 if (mod(nn-1,2).eq.0) then
1726 j = ((nn-1)/2)+1
1727
1728 ! Odd spectral polynomial degree
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(k) &
1736 + beta13(ie)*ct(j) +gamma1(ie)*ct(j)*ct(k)
1737 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1738 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1739 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1740 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1741
1742 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1743 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1744 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1745 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1746 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1747 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1748
1749 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1750 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1751 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1752 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1753 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1754 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
1755
1756 det_j = dxdz * (dydx*dzdy - dzdx*dydy) &
1757 - dydz * (dxdx*dzdy - dzdx*dxdy) &
1758 + dzdz * (dxdx*dydy - dydx*dxdy)
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) * ww(i) * ww(k) * dabs(det_j) / ellez
1764 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in -1) +3)) + term
1765 enddo
1766 enddo
1767
1768
1769
1770 !****** Situation S5 or S6 - begin ******
1771
1772 elseif (sum_node_bottom.eq.18) then
1773
1774 !*** Situation S5 ***
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 !*** Situation S6 ***
1783
1784 else
1785
1786 sit = 2
1787
1788 endif
1789
1790 ! Even spectral polynomial degree
1791 if (mod(nn-1,2).eq.0) then
1792 i = ((nn-1)/2)+1
1793
1794 ! Odd spectral polynomial degree
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)*ct(j)*ct(k)
1804 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1805 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1806 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1807 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1808
1809 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1810 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1811 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1812 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1813 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1814 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1815
1816 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1817 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1818 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1819 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1820 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1821 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
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) * ww(j) * ww(k) * dabs(det_j) / ellez
1831 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in -1) +3)) + term
1832 enddo
1833 enddo
1834
1835 endif
1836
1837
1838 endif
1839 endif !Check on the Plane Wave Material - end
1840
1841 enddo
1842 endif ! Plane Wave Z - end
1843
1844
1845 !############################################################################
1846 !################## PLANE WAVE END ##############
1847 !############################################################################
1848
1849
1850 !############################################################################
1851 !################## SEISMIC MOMENT BEGIN ##############
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(j)*ct(k)
1872 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1873 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1874 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1875 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1876
1877 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1878 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1879 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1880 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1881 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1882 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1883
1884 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1885 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1886 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1887 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1888 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1889 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
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(ipsism,isism)) then
1900 facsmom(isism,1) = facsmom(isism,1) + det_j * ww(i) * ww(j) * ww(k)
1901 facsmom(isism,2) = facsmom(isism,2) + det_j * ww(i) * ww(j) * ww(k)
1902 facsmom(isism,3) = facsmom(isism,3) + det_j * ww(i) * ww(j) * ww(k)
1903 facsmom(isism,4) = facsmom(isism,4) + det_j * ww(i) * ww(j) * ww(k)
1904 facsmom(isism,5) = facsmom(isism,5) + det_j * ww(i) * ww(j) * ww(k)
1905 facsmom(isism,6) = facsmom(isism,6) + det_j * ww(i) * ww(j) * ww(k)
1906
1907 length_cns = length_cns + 1
1908
1909 endif
1910 enddo !i
1911 enddo !j
1912 enddo !k
1913 endif !if (fn.gt.0) then
1914
1915 enddo !ip
1916
1917 enddo !isism
1918
1919 endif !if (nl_sism.gt.0) then
1920
1921
1922 !############################################################################
1923 !################## SEISMIC MOMENT END ##############
1924 !############################################################################
1925
1926
1927 !############################################################################
1928 !################## EXPLOSIVE SOURCE BEGIN ##############
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(j)*ct(k)
1946 dydx = alfa21(ie) +beta22(ie)*ct(k) &
1947 + beta23(ie)*ct(j) +gamma2(ie)*ct(j)*ct(k)
1948 dzdx = alfa31(ie) +beta32(ie)*ct(k) &
1949 + beta33(ie)*ct(j) +gamma3(ie)*ct(j)*ct(k)
1950
1951 dxdy = alfa12(ie) +beta11(ie)*ct(k) &
1952 + beta13(ie)*ct(i) +gamma1(ie)*ct(k)*ct(i)
1953 dydy = alfa22(ie) +beta21(ie)*ct(k) &
1954 + beta23(ie)*ct(i) +gamma2(ie)*ct(k)*ct(i)
1955 dzdy = alfa32(ie) +beta31(ie)*ct(k) &
1956 + beta33(ie)*ct(i) +gamma3(ie)*ct(k)*ct(i)
1957
1958 dxdz = alfa13(ie) +beta11(ie)*ct(j) &
1959 + beta12(ie)*ct(i) +gamma1(ie)*ct(i)*ct(j)
1960 dydz = alfa23(ie) +beta21(ie)*ct(j) &
1961 + beta22(ie)*ct(i) +gamma2(ie)*ct(i)*ct(j)
1962 dzdz = alfa33(ie) +beta31(ie)*ct(j) &
1963 + beta32(ie)*ct(i) +gamma3(ie)*ct(i)*ct(j)
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(ipexpl,iexpl)) then
1973 facsexpl(iexpl,1) = facsexpl(iexpl,1) + det_j * ww(i) * ww(j) * ww(k)
1974 facsexpl(iexpl,2) = facsexpl(iexpl,2) + det_j * ww(i) * ww(j) * ww(k)
1975 facsexpl(iexpl,3) = facsexpl(iexpl,3) + det_j * ww(i) * ww(j) * ww(k)
1976 facsexpl(iexpl,4) = facsexpl(iexpl,4) + det_j * ww(i) * ww(j) * ww(k)
1977 facsexpl(iexpl,5) = facsexpl(iexpl,5) + det_j * ww(i) * ww(j) * ww(k)
1978 facsexpl(iexpl,6) = facsexpl(iexpl,6) + det_j * ww(i) * ww(j) * ww(k)
1979
1980 length_cne = length_cne + 1
1981
1982 endif
1983 enddo !i
1984 enddo !j
1985 enddo !k
1986 endif !if (fn.gt.0) then
1987
1988 enddo !ip
1989
1990 enddo !iexpl
1991
1992
1993 endif !if (nl_expl.gt.0) then
1994
1995
1996
1997 !############################################################################
1998 !################## EXPLOSIVE SOURCE END ##############
1999 !############################################################################
2000 deallocate(ct,dd,ww)
2001 enddo
2002
2003
2004
2005
2006!############################################################################
2007!################## SEISMIC MOMENT BEGIN ##############
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, mpi_sum, mpi_comm, mpi_ierr)
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 ! Seismic Tensor is given is terms of norma and slip vector:
2038 ! nx, ny, nz and sx, sy, sz
2039 ! i.e.: SULMONA Case
2040 if(facsmom(isism,1) .ne. 0) &
2041 facsmom(isism,1) = 1/facsmom(isism,1) &
2042 * (slip1_sism*norm1_sism+slip1_sism*norm1_sism) &
2043 * amp_sism
2044 if(facsmom(isism,2) .ne. 0) &
2045 facsmom(isism,2) = 1/facsmom(isism,2) &
2046 * (slip2_sism*norm2_sism+slip2_sism*norm2_sism) &
2047 * amp_sism
2048 if(facsmom(isism,3) .ne. 0) &
2049 facsmom(isism,3) = 1/facsmom(isism,3) &
2050 * (slip3_sism*norm3_sism+slip3_sism*norm3_sism) &
2051 * amp_sism
2052 if(facsmom(isism,4) .ne. 0) &
2053 facsmom(isism,4) = 1/facsmom(isism,4) &
2054 * (slip2_sism*norm3_sism+slip3_sism*norm2_sism) &
2055 * amp_sism
2056 if(facsmom(isism,5) .ne. 0) &
2057 facsmom(isism,5) = 1/facsmom(isism,5) &
2058 * (slip1_sism*norm3_sism+slip3_sism*norm1_sism) &
2059 * amp_sism
2060 if(facsmom(isism,6) .ne. 0) &
2061 facsmom(isism,6) = 1/facsmom(isism,6) &
2062 * (slip1_sism*norm2_sism+slip2_sism*norm1_sism) &
2063 * amp_sism
2064
2065 tausmom(isism,1) = tau_sism
2066
2067 ! Seismic Tensor is given is terms of norma and slip vector
2068 ! mxx,myy,mzz,mxy,mxz,myz
2069 ! i.e.: CASHIMA (ONLY)
2070 !facsmom(isism,1) = 1/facsmom(isism,1) &
2071 ! * slip1_sism
2072 !facsmom(isism,2) = 1/facsmom(isism,2) &
2073 ! * slip2_sism
2074 !facsmom(isism,3) = 1/facsmom(isism,3) &
2075 ! * slip3_sism
2076 !facsmom(isism,4) = 1/facsmom(isism,4) &
2077 ! * norm1_sism
2078 !facsmom(isism,5) = 1/facsmom(isism,5) &
2079 ! * norm2_sism
2080 !facsmom(isism,6) = 1/facsmom(isism,6) &
2081 ! * norm3_sism
2082 !facsmom(isism,1) = 1/facsmom(isism,1) &
2083 ! * 1 &
2084 ! * amp_sism
2085 !facsmom(isism,2) = 1/facsmom(isism,2) &
2086 ! * 1 &
2087 ! * amp_sism
2088 !facsmom(isism,3) = 1/facsmom(isism,3) &
2089 ! * 1 &
2090 ! * amp_sism
2091 !facsmom(isism,4) = 0.0d0
2092 !facsmom(isism,5) = 0.0d0
2093 !facsmom(isism,6) = 0.0d0
2094 enddo
2095
2096 endif
2097
2098
2099!############################################################################
2100!################## SEISMIC MOMENT END ##############
2101!############################################################################
2102
2103
2104!############################################################################
2105!################## EXPLOSIVE SOURCE BEGIN ##############
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, mpi_sum, mpi_comm, mpi_ierr)
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!################## EXPLOSIVE SOURCE END ##############
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))
2164 call make_lgl_nw(nn,ct,ww,dd)
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 ! Test load
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(im,3)
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 !EXTERNAL LOAD PAPER WITH BLANCA
2218 fmat(fn,(3*(in -1) +1)) = fmat(fn,(3*(in -1) +1)) + det_j * ww(i) * ww(j) * ww(k) * ( &
2219 4.d0*pi**2.d0*dcos(pi*y)*dcos(pi*z)*dsin(pi*y)*dsin(pi*z) * &
2220 (2.d0*lambda - 8.d0*mu + 9.d0*rho - 4.d0*lambda*dcos(pi*x)**2.d0 &
2221 + 8.d0*mu*dcos(pi*x)**2.d0 - 9.d0*rho*dcos(pi*x)**2.d0) )
2222
2223 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1) +2)) + det_j * ww(i) * ww(j) * ww(k) * ( &
2224 4.d0*pi**2.d0*dcos(pi*x)*dcos(pi*z)*dsin(pi*x)*dsin(pi*z) * &
2225 (2.d0*lambda + 12.d0*mu - 9.d0*rho - 4.d0*lambda*dcos(pi*y)**2.d0 &
2226 - 16.d0*mu*dcos(pi*y)**2.d0 + 9.d0*rho*dcos(pi*y)**2.d0) )
2227
2228 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in -1) +3)) + det_j * ww(i) * ww(j) * ww(k) * (&
2229 4.d0*pi**2.d0*dcos(pi*x)*dcos(pi*y)*dsin(pi*x)*dsin(pi*y) * &
2230 (2.d0*lambda + 12.d0*mu - 9.d0*rho - 4.d0*lambda*dcos(pi*z)**2.d0 &
2231 - 16.d0*mu*dcos(pi*z)**2.d0 + 9.d0*rho*dcos(pi*z)**2.d0) )
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!if ((cs_bc_loc(ie) - cs_bc_loc(ie -1) -1) .ne. nn*nn) then
2249!deallocate(ct,ww,dd)
2250nn = int(sqrt(real(cs_bc_loc(ie) - cs_bc_loc(ie -1) -1)))
2251
2252allocate(ct(nn),ww(nn),dd(nn,nn))
2253call MAKE_LGL_NW(nn,ct,ww,dd)
2254!endif
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 ! Neumann load X
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 ! Neumann load Y
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 ! Neumann load Z
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; !val_neuN(il,1)
2412 v2 = 0 !val_neuN(il,2)
2413 v3 = 0 !val_neuN(il,3)
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 !write(*,*) area
2432 !read(*,*)
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(k1*x + k2*y + k3*z)
2436 v2 = k3*lambda*n2*cos(k1*x + k2*y + k3*z) + k2*mu*n3*cos(k1*x + k2*y + k3*z)
2437 v3 = n3*(k3*lambda*cos(k1*x + k2*y + k3*z) + 2*k3*mu*cos(k1*x + k2*y + k3*z)) &
2438 + k1*mu*n1*cos(k1*x + k2*y + k3*z) + k2*mu*n2*cos(k1*x + k2*y + k3*z)
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(k1*x + k2*y + k3*z)
2444 v2 = - k3*lambda*n2*sin(k1*x + k2*y + k3*z) - k2*mu*n3*sin(k1*x + k2*y + k3*z)
2445 v3 = - n3*(k3*lambda*sin(k1*x + k2*y + k3*z) + 2*k3*mu*sin(k1*x + k2*y + k3*z)) &
2446 - k1*mu*n1*sin(k1*x + k2*y + k3*z) - k2*mu*n2*sin(k1*x + k2*y + k3*z)
2447
2448 v1 = -v1; v2 = -v2; v3 = -v3;
2449 else
2450 write(*,*) 'Case not implemented!!! Check inputs.'
2451
2452 endif
2453 ! v1 = 0; v2 = 0; v3 = 0;
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(2*pi*z) &
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*pi*z) &
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*pi*z) &
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.25*area * v1 * ww(i)*ww(j)
2483 fmat(fn,(3*(in -1) +2)) = fmat(fn,(3*(in -1) +2)) + 0.25*area * v2 * ww(i)*ww(j)
2484 fmat(fn,(3*(in -1) +3)) = fmat(fn,(3*(in -1) +3)) + 0.25*area * v3 * ww(i)*ww(j)
2485
2486 enddo
2487 enddo
2488 endif
2489 endif
2490enddo
2491
2492endif
2493
2494! ------------------------------------------------------------------------
2495
2496if (nl_dirx.gt.0) then ! Dirichlet X
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 ! Dirichlet Y
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 ! Dirichlet Z
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, normal_nz_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! deallocate(ct,ww,dd)
2599
2600! do i = 1, nfunc
2601! do j = 1, nnod_loc
2602! write(*,*) fmat(i,3*(j-1)+1), fmat(i,3*(j-1)+2),fmat(i,3*(j-1)+3)
2603! enddo
2604! read(*,*)
2605! enddo
2606!stop
2607
2608
2609return
2610end subroutine make_extint_forces
2611
2612
2613
get_elem_from_surf
subroutine get_elem_from_surf(nnz_bc, cs_bc, n1, n2, n3, n4, ie_surf)
Find element index from verteces number.
Definition GET_ELEM_FROM_SURF.f90:33
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_mech_prop_nh_enhanced
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
Definition GET_MECH_PROP_NH_ENHANCED.f90:28
get_minvalues
subroutine get_minvalues(i_glo, v_glo, n_glo, i_loc, n_loc, np)
Computes positions of minimum values.
Definition GET_MINVALUES.f90:32
get_nearest_node
subroutine get_nearest_node(n, xs, ys, zs, xt, yt, zt, nt, dist_min)
Computes the nearest node with respet to (xt,yt,zt).
Definition GET_NEAREST_NODE.f90:35
get_node_from_face
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.
Definition GET_NODE_FROM_FACE.f90:34
make_extint_forces
subroutine make_extint_forces(nnod_loc, xs_loc, ys_loc, zs_loc, local
Computes external loads.
Definition MAKE_EXTINT_FORCES.f90:96
make_lgl_nw
subroutine make_lgl_nw(nb_pnt, xq, wq, dd)
Makes Gauss-Legendre-Lobatto nodes, weigths and spectral derivatives.
Definition MAKE_LGL_NW.f90:29
make_normal
subroutine make_normal(ind, xs1, xs2, xs3, xs4, ys1, ys2, ys3, ys4, zs1, zs2, zs3, zs4, nx, ny, nz, par, arr)
Makes normal vector of a given surface.
Definition MAKE_NORMAL.f90:48
phi
real *8 function phi(ng, csii, csij)
Evaluates Lengendre basis function on a specific point.
Definition PHI.f90:31
speed_exit_codes
SPEED exit codes.
Definition MODULES.f90:25
speed_exit_codes::exit_surf_notfound
integer, parameter exit_surf_notfound
Definition MODULES.f90:35