All Classes Files Functions Variables Pages
module_receiver.f90
Go to the documentation of this file.
1 !=====================================================================================================================================
2 ! EFISPEC3D !
3 ! (Elements FInis SPECtraux 3D) !
4 ! !
5 ! http://efispec.free.fr !
6 ! !
7 ! !
8 ! This file is part of EFISPEC3D !
9 ! Please refer to http://efispec.free.fr if you use it or part of it !
10 ! !
11 ! !
12 !1 ---> French License: CeCILL V2 !
13 ! !
14 ! Copyright BRGM 2009 contributeurs : Florent DE MARTIN !
15 ! David MICHEA !
16 ! Philippe THIERRY !
17 ! !
18 ! Contact: f.demartin at brgm.fr !
19 ! !
20 ! Ce logiciel est un programme informatique servant a resoudre l'equation du !
21 ! mouvement en trois dimensions via une methode des elements finis spectraux. !
22 ! !
23 ! Ce logiciel est regi par la licence CeCILL soumise au droit francais et !
24 ! respectant les principes de diffusion des logiciels libres. Vous pouvez !
25 ! utiliser, modifier et/ou redistribuer ce programme sous les conditions de la !
26 ! licence CeCILL telle que diffusee par le CEA, le CNRS et l'INRIA sur le site !
27 ! "http://www.cecill.info". !
28 ! !
29 ! En contrepartie de l'accessibilite au code source et des droits de copie, de !
30 ! modification et de redistribution accordes par cette licence, il n'est offert !
31 ! aux utilisateurs qu'une garantie limitee. Pour les memes raisons, seule une !
32 ! responsabilite restreinte pese sur l'auteur du programme, le titulaire des !
33 ! droits patrimoniaux et les concedants successifs. !
34 ! !
35 ! A cet egard l'attention de l'utilisateur est attiree sur les risques associes !
36 ! au chargement, a l'utilisation, a la modification et/ou au developpement et a !
37 ! la reproduction du logiciel par l'utilisateur etant donne sa specificite de !
38 ! logiciel libre, qui peut le rendre complexe a manipuler et qui le reserve donc !
39 ! a des developpeurs et des professionnels avertis possedant des connaissances !
40 ! informatiques approfondies. Les utilisateurs sont donc invites a charger et !
41 ! tester l'adequation du logiciel a leurs besoins dans des conditions permettant !
42 ! d'assurer la securite de leurs systemes et ou de leurs donnees et, plus !
43 ! generalement, a l'utiliser et l'exploiter dans les memes conditions de !
44 ! securite. !
45 ! !
46 ! Le fait que vous puissiez acceder a cet en-tete signifie que vous avez pris !
47 ! connaissance de la licence CeCILL et que vous en avez accepte les termes. !
48 ! !
49 ! !
50 !2 ---> International license: GNU GPL V3 !
51 ! !
52 ! EFISPEC3D is a computer program that solves the three-dimensional equations of !
53 ! motion using a finite spectral-element method. !
54 ! !
55 ! Copyright (C) 2009 Florent DE MARTIN !
56 ! !
57 ! Contact: f.demartin at brgm.fr !
58 ! !
59 ! This program is free software: you can redistribute it and/or modify it under !
60 ! the terms of the GNU General Public License as published by the Free Software !
61 ! Foundation, either version 3 of the License, or (at your option) any later !
62 ! version. !
63 ! !
64 ! This program is distributed in the hope that it will be useful, but WITHOUT ANY !
65 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A !
66 ! PARTICULAR PURPOSE. See the GNU General Public License for more details. !
67 ! !
68 ! You should have received a copy of the GNU General Public License along with !
69 ! this program. If not, see http://www.gnu.org/licenses/. !
70 ! !
71 ! !
72 !3 ---> Third party libraries !
73 ! !
74 ! EFISPEC3D uses the following of third party libraries: !
75 ! !
76 ! --> METIS 5.1.0 !
77 ! see http://glaros.dtc.umn.edu/gkhome/metis/metis/overview !
78 ! !
79 ! --> Lib_VTK_IO !
80 ! see S. Zaghi's website: https://github.com/szaghi/Lib_VTK_IO !
81 ! !
82 ! --> INTERP_LINEAR !
83 ! see J. Burkardt website: http://people.sc.fsu.edu/~jburkardt/ !
84 ! !
85 ! --> SAC !
86 ! http://ds.iris.edu/files/sac-manual/ !
87 ! !
88 ! --> EXODUS II !
89 ! http://sourceforge.net/projects/exodusii/ !
90 ! !
91 ! --> NETCDF !
92 ! http://www.unidata.ucar.edu/software/netcdf/ !
93 ! !
94 ! --> HDF5 !
95 ! http://www.hdfgroup.org/HDF5/ !
96 ! !
97 ! Some of these libraries are located in directory lib and pre-compiled !
98 ! with intel compiler for x86_64 platform. !
99 ! !
100 ! !
101 !4 ---> Related Articles !
102 ! !
103 ! De Martin, F., Matsushima, M., Kawase, H. (BSSA, 2013) !
104 ! Impact of geometric effects on near-surface Green's functions !
105 ! doi:10.1785/0120130039 !
106 ! !
107 ! Aochi, H., Ducellier, A., Dupros, F., Delatre, M., Ulrich, T., De Martin, F., !
108 ! Yoshimi, M., (Pure Appl. Geophys. 2013) !
109 ! Finite Difference Simulations of Seismic Wave Propagation for the 2007 Mw 6.6 !
110 ! Niigata-ken Chuetsu-Oki Earthquake: Validity of Models and Reliable !
111 ! Input Ground Motion in the Near-Field !
112 ! doi:10.1007/s00024-011-0429-5 !
113 ! !
114 ! De Martin, F. (BSSA, 2011) !
115 ! Verification of a Spectral-Element Method Code for the Southern California !
116 ! Earthquake Center LOH.3 Viscoelastic Case !
117 ! doi:10.1785/0120100305 !
118 ! !
119 !=====================================================================================================================================
120 
123 
126 module mod_receiver
127 
128  use mpi
129 
130  implicit none
131 
132  private
133 
134  public :: init_hexa_receiver
135  public :: init_quad_receiver
136  public :: search_closest_hexa_gll
137  public :: search_closest_quad_gll
138  public :: compute_info_hexa_receiver
139  public :: compute_info_quad_receiver
140  public :: write_receiver_output
141  public :: is_inside_quad
142  private :: is_inside_polygon
143 
144  contains
145 
146 !
147 !
153 !***********************************************************************************************************************************************************************************
154  subroutine init_hexa_receiver()
155 !***********************************************************************************************************************************************************************************
156 
157  use mpi
158 
159  use mod_global_variables, only :&
160  cg_prefix&
161  ,tg_receiver_hexa&
162  ,ig_nreceiver_hexa&
163  ,ig_lst_unit&
164  ,ig_myrank&
165  ,ig_ncpu&
166  ,type_receiver_hexa&
167  ,ig_receiver_saving_incr&
168  ,lg_output_debug_file&
169  ,cg_myrank&
170  ,ig_hexa_receiver_unit&
171  ,get_newunit&
172  ,ig_ndt
173 
174  implicit none
175 
176  real , dimension(ig_ncpu) :: rldum1
177  real , dimension(ig_ncpu) :: maxdmin_all_cpu
178  real :: maxdmin
179  real :: rldmin
180 
181  integer, dimension(1) :: min_loc
182  integer, dimension(ig_ncpu) :: ilrcpu
183  integer :: ios
184  integer :: irec
185  integer :: j
186  integer :: ilnrec
187 
188  character(len=1) :: ctmp
189  character(len=6) :: crec
190 
191  type(type_receiver_hexa), allocatable, dimension(:) :: tlrinf
192 
193 !
194 !---->initialize total number of receivers in cpu myrank
195  ig_nreceiver_hexa = 0
196 !
197 !
198 !---->if file *.vor exists, then we read receiver x,y,z location
199  open(unit=100,file=trim(cg_prefix)//".vor",status='old',iostat=ios)
200 
201  if( (ios == 0) .and. (ig_receiver_saving_incr <= ig_ndt) ) then
202 
203  if (ig_myrank == 0) then
204  write(ig_lst_unit,'(a)') " "
205  write(ig_lst_unit,'(a)') "information about receivers in file "//trim(cg_prefix)//".vor"
206  endif
207 
208  do while (.true.)
209  read(unit=100,fmt=*,iostat=ios) ctmp
210  if (ios /= 0) exit
211  ig_nreceiver_hexa = ig_nreceiver_hexa + 1
212  enddo
213  rewind(100)
214 
215  allocate(tlrinf(ig_nreceiver_hexa))
216  ilnrec = 0
217 
218  do irec = 1,ig_nreceiver_hexa
219 
220 !
221 !---------->initialize tlrinf
222  tlrinf(irec)%x = 0.0
223  tlrinf(irec)%y = 0.0
224  tlrinf(irec)%z = 0.0
225  tlrinf(irec)%xi = 0.0
226  tlrinf(irec)%eta = 0.0
227  tlrinf(irec)%zeta = 0.0
228  tlrinf(irec)%dmin = huge(rldmin)
229  tlrinf(irec)%lag(:,:,:) = 0.0
230  tlrinf(irec)%gll(:,:,:) = 0
231  tlrinf(irec)%cpu = 0
232  tlrinf(irec)%iel = 0
233  tlrinf(irec)%kgll = 0
234  tlrinf(irec)%lgll = 0
235  tlrinf(irec)%mgll = 0
236  tlrinf(irec)%rglo = 0
237 
238  read(unit=100,fmt=*) tlrinf(irec)%x,tlrinf(irec)%y,tlrinf(irec)%z
239 !
240 !---------->find the closest element and gll point to the receiver irec
241  call search_closest_hexa_gll(tlrinf(irec)%x ,tlrinf(irec)%y ,tlrinf(irec)%z ,tlrinf(irec)%dmin&
242  ,tlrinf(irec)%kgll,tlrinf(irec)%lgll,tlrinf(irec)%mgll,tlrinf(irec)%iel )
243 !
244 !---------->check to which cpu the receiver belongs (i.e., cpu which has the smallest dmin)
245  call mpi_allgather(tlrinf(irec)%dmin,1,mpi_real,rldum1,1,mpi_real,mpi_comm_world,ios)
246  min_loc = minloc(rldum1(1:ig_ncpu))
247 
248  if (min_loc(1) == ig_myrank+1) then
249  tlrinf(irec)%cpu = ig_myrank + 1
250  ilnrec = ilnrec + 1
251  else
252  tlrinf(irec)%cpu = min_loc(1)
253  endif
254 
255  enddo
256  close(100)
257 !
258 !------->transfer array tlrinf which knows all the receivers to array tg_receiver_hexa that needs only to know the receiver of cpu myrank
259  if (ilnrec > 0) then
260  allocate(tg_receiver_hexa(ilnrec))
261  j = 0
262  do irec = 1,ig_nreceiver_hexa
263  if (tlrinf(irec)%cpu == ig_myrank+1) then
264  j = j + 1
265  tg_receiver_hexa(j) = tlrinf(irec)
266  tg_receiver_hexa(j)%rglo = irec
267  endif
268  enddo
269  ig_nreceiver_hexa = ilnrec
270  else
271  ig_nreceiver_hexa = 0
272  endif
273  deallocate(tlrinf)
274 !
275 !------->cpu 0 gather the number of receiver (ig_nreceiver_hexa) of the other cpus
276  call mpi_gather(ig_nreceiver_hexa,1,mpi_integer,ilrcpu,1,mpi_integer,0,mpi_comm_world,ios)
277 !
278 !------->compute information for receivers that belongs to cpu myrank and send 'dmin' to cpu 0
279  do irec = 1,ig_nreceiver_hexa
280 
281  call compute_info_hexa_receiver(tg_receiver_hexa(irec))
282 
283  enddo
284 
285 !
286 !------->all cpus compute their maximum dmin
287  if (ig_nreceiver_hexa > 0) then
288 
289  maxdmin = -huge(maxdmin)
290 
291  do irec = 1,ig_nreceiver_hexa
292  maxdmin = max(maxdmin,tg_receiver_hexa(irec)%dmin)
293  enddo
294 
295  else
296 
297  maxdmin = 0.0
298 
299  endif
300 
301 !
302 !------->cpu 0 gathers dmin of all cpus and write the maximum dmin in *.lst
303  call mpi_gather(maxdmin,1,mpi_real,maxdmin_all_cpu,1,mpi_real,0,mpi_comm_world,ios)
304 
305  if (ig_myrank == 0) then
306 
307  write(ig_lst_unit,'(a,e14.7)') "maximum localisation error of all receivers in free surface = ",maxval(maxdmin_all_cpu)
308  call flush(ig_lst_unit)
309 
310  endif
311 
312 !
313 !------->open receivers file once and for all
314  if ( ig_nreceiver_hexa > 0 ) then
315 
316  allocate(ig_hexa_receiver_unit(ig_nreceiver_hexa))
317 
318  do irec = 1,ig_nreceiver_hexa
319 
320  write(crec,'(i6.6)') tg_receiver_hexa(irec)%rglo
321 
322  open(unit = get_newunit(ig_hexa_receiver_unit(irec))&
323  ,file = trim(cg_prefix)//".vor."//trim(adjustl(crec))//".gpl"&
324  ,status = 'replace'&
325  ,action = 'write'&
326  ,access = 'sequential'&
327  ,form = 'unformatted'&
328  ,buffered = 'yes'&
329  ,recordtype = 'stream')
330 
331  enddo
332 
333  endif
334 
335  if (lg_output_debug_file) then
336  open(unit=100,file="debug.cpu."//trim(cg_myrank)//".volume.receivers.info")
337  do irec = 1,ig_nreceiver_hexa
338  write(unit=100,fmt='(/,a,i10 )') "Receiver ",tg_receiver_hexa(irec)%rglo
339  write(unit=100,fmt='( a,i10 )') " in core ",tg_receiver_hexa(irec)%cpu-1
340  write(unit=100,fmt='( a,i10 )') " in hexa ",tg_receiver_hexa(irec)%iel
341  write(unit=100,fmt='( 3a )') " X "," Y "," Z "
342  write(unit=100,fmt='(3(e14.7,1x))') tg_receiver_hexa(irec)%x,tg_receiver_hexa(irec)%y,tg_receiver_hexa(irec)%z
343  write(unit=100,fmt='( 3a )') " XI "," ETA "," ZETA "
344  write(unit=100,fmt='(3(e14.7,1x))') tg_receiver_hexa(irec)%xi,tg_receiver_hexa(irec)%eta,tg_receiver_hexa(irec)%zeta
345  enddo
346  close(100)
347  endif
348 
349  else
350 
351  if (ig_myrank == 0) then
352  write(ig_lst_unit,'(/a)') "no volume receiver computed"
353  call flush(ig_lst_unit)
354  endif
355 
356  endif
357 
358  return
359 !***********************************************************************************************************************************************************************************
360  end subroutine init_hexa_receiver
361 !***********************************************************************************************************************************************************************************
362 
363 !
364 !
370 !***********************************************************************************************************************************************************************************
371  subroutine init_quad_receiver()
372 !***********************************************************************************************************************************************************************************
373 
374  use mpi
375 
376  use mod_global_variables, only :&
377  cg_prefix&
378  ,error_stop&
379  ,tg_receiver_quad&
380  ,ig_nquad_fsurf&
381  ,ig_quad_nnode&
382  ,ig_quadf_gll_glonum&
383  ,ig_quadf_gnode_glonum&
384  ,ig_nreceiver_quad&
385  ,ig_lst_unit&
386  ,ig_myrank&
387  ,ig_ncpu&
388  ,type_receiver_quad&
389  ,ig_receiver_saving_incr&
390  ,lg_output_debug_file&
391  ,cg_myrank&
392  ,ig_quad_receiver_unit&
393  ,get_newunit&
394  ,ig_ndt
395 
396  implicit none
397 
398  real , dimension(ig_ncpu) :: maxdmin_all_cpu
399  real :: maxdmin
400  real :: rldmin
401 
402  integer, dimension(ig_ncpu) :: ilrcpu
403  integer :: ios
404  integer :: icpu
405  integer :: iloc
406  integer :: irec
407  integer :: j
408  integer :: ilnrec
409 
410  character(len= 1) :: ctmp
411  character(len= 6) :: crec
412  character(len=255) :: info
413 
414  logical :: is_inside
415  logical, dimension(ig_ncpu) :: is_inside_all_cpu
416 
417  type(type_receiver_quad), allocatable, dimension(:) :: tlrinf
418 
419 !
420 !---->initialize total number of receivers in cpu myrank
421  ig_nreceiver_quad = 0
422 !
423 !
424 !---->if file *receiver exists, then we read receiver x,y,z location
425  open(unit=100,file=trim(cg_prefix)//".fsr",status='old',iostat=ios)
426 
427  if( (ios == 0) .and. (ig_receiver_saving_incr <= ig_ndt) ) then
428 
429  if (ig_myrank == 0) then
430  write(ig_lst_unit,'(a)') " "
431  write(ig_lst_unit,'(a)') "information about receivers in file "//trim(cg_prefix)//".fsr"
432  endif
433 
434  do while (.true.)
435  read(unit=100,fmt=*,iostat=ios) ctmp
436  if (ios /= 0) exit
437  ig_nreceiver_quad = ig_nreceiver_quad + 1
438  enddo
439  rewind(100)
440 
441  allocate(tlrinf(ig_nreceiver_quad))
442  ilnrec = 0
443 
444  do irec = 1,ig_nreceiver_quad
445 
446 !
447 !---------->initialize tlrinf
448  tlrinf(irec)%x = 0.0
449  tlrinf(irec)%y = 0.0
450  tlrinf(irec)%z = 0.0
451  tlrinf(irec)%xi = 0.0
452  tlrinf(irec)%eta = 0.0
453  tlrinf(irec)%dmin = huge(rldmin)
454  tlrinf(irec)%lag(:,:) = 0.0
455  tlrinf(irec)%gll(:,:) = 0
456  tlrinf(irec)%cpu = -1
457  tlrinf(irec)%iel = 0
458  tlrinf(irec)%lgll = 0
459  tlrinf(irec)%mgll = 0
460  tlrinf(irec)%rglo = 0
461 
462  read(unit=100,fmt=*) tlrinf(irec)%x,tlrinf(irec)%y
463 
464 !
465 !---------->find which quadrangle element contains the receiver irec
466  if (ig_nquad_fsurf > 0) then
467 
468  call is_inside_quad(tlrinf(irec)%x,tlrinf(irec)%y,ig_quadf_gnode_glonum,ig_nquad_fsurf,ig_quad_nnode,tlrinf(irec)%dmin,tlrinf(irec)%lgll,tlrinf(irec)%mgll,tlrinf(irec)%iel,is_inside)
469 
470  else
471 
472  is_inside = .false.
473 
474  endif
475 
476 !
477 !---------->cpu0 gather all is_inside to check: 1) if the receiver is inside a quadrangle element: if not abort computation 2) if the receiver is found by multiple cpus: if yes, affect the receiver to only one cpu
478  call mpi_gather(is_inside,1,mpi_logical,is_inside_all_cpu,1,mpi_logical,0,mpi_comm_world,ios)
479 
480  if (ig_myrank == 0) then
481 !
482 !------------->search the first cpu that contains the receiver
483  iloc = -1
484  do icpu = 1,ig_ncpu
485 
486  if (is_inside_all_cpu(icpu)) then
487  iloc = icpu
488  exit
489  endif
490 
491  enddo
492 
493 !
494 !------------->if the receiver is not inside any quadrangle element among all the cpus --> abort computation
495  if (iloc == -1) then
496 
497  write(info,'(a)') "error in subroutine init_quad_receiver: receiver not found over all cpus"
498  call error_stop(info)
499 
500  else
501 
502  do icpu = iloc+1,ig_ncpu
503  is_inside_all_cpu(icpu) = .false.
504  enddo
505 
506  endif
507 
508  endif
509 
510 !
511 !---------->cpu0 scatter array is_inside_all_cpu. The cpu that obtain the value is_inside = .true. will compute the receiver irec
512  call mpi_scatter(is_inside_all_cpu,1,mpi_logical,is_inside,1,mpi_logical,0,mpi_comm_world,ios)
513 
514  if (is_inside) then
515 
516  tlrinf(irec)%cpu = ig_myrank + 1
517  ilnrec = ilnrec + 1
518 
519  endif
520 
521  enddo
522 
523  close(100)
524 
525 !
526 !------->transfer array tlrinf which knows all the receivers to array tg_receiver_quad that needs only to know the receiver of cpu myrank
527  if (ilnrec > 0) then
528 
529  allocate(tg_receiver_quad(ilnrec))
530 
531  j = 0
532  do irec = 1,ig_nreceiver_quad
533  if (tlrinf(irec)%cpu == ig_myrank+1) then
534  j = j + 1
535  tg_receiver_quad(j) = tlrinf(irec)
536  tg_receiver_quad(j)%rglo = irec
537  endif
538  enddo
539 
540  ig_nreceiver_quad = ilnrec
541 
542  else
543 
544  ig_nreceiver_quad = 0
545 
546  endif
547 
548  deallocate(tlrinf)
549 !
550 !------->cpu 0 gather the number of receiver (ig_nreceiver_quad) of the other cpus
551  call mpi_gather(ig_nreceiver_quad,1,mpi_integer,ilrcpu,1,mpi_integer,0,mpi_comm_world,ios)
552 !
553 !------->compute information for receivers that belongs to cpu myrank and send 'dmin' to cpu 0
554  do irec = 1,ig_nreceiver_quad
555 
556  call compute_info_quad_receiver(tg_receiver_quad(irec))
557 
558  enddo
559 
560 !
561 !------->all cpus compute their maximum dmin
562  if (ig_nreceiver_quad > 0) then
563 
564  maxdmin = -huge(maxdmin)
565 
566  do irec = 1,ig_nreceiver_quad
567  maxdmin = max(maxdmin,tg_receiver_quad(irec)%dmin)
568  enddo
569 
570  else
571 
572  maxdmin = 0.0
573 
574  endif
575 
576 !
577 !------->cpu 0 gathers dmin of all cpus and write the maximum dmin in *.lst
578  call mpi_gather(maxdmin,1,mpi_real,maxdmin_all_cpu,1,mpi_real,0,mpi_comm_world,ios)
579 
580  if (ig_myrank == 0) then
581 
582  write(ig_lst_unit,'(a,e14.7)') "maximum localisation error of all receivers in free surface = ",maxval(maxdmin_all_cpu)
583  call flush(ig_lst_unit)
584 
585  endif
586 !
587 !------->open receivers file once and for all
588  if ( ig_nreceiver_quad > 0 ) then
589 
590  allocate(ig_quad_receiver_unit(ig_nreceiver_quad))
591 
592  do irec = 1,ig_nreceiver_quad
593 
594  write(crec,'(i6.6)') tg_receiver_quad(irec)%rglo
595 
596  open(unit = get_newunit(ig_quad_receiver_unit(irec))&
597  ,file = trim(cg_prefix)//".fsr."//trim(adjustl(crec))//".gpl"&
598  ,status = 'replace'&
599  ,action = 'write'&
600  ,access = 'sequential'&
601  ,form = 'unformatted'&
602  ,buffered = 'yes'&
603  ,recordtype = 'stream')
604 
605  enddo
606 
607  endif
608 
609  if (lg_output_debug_file) then
610  open(unit=100,file="debug.cpu."//trim(cg_myrank)//".freesurface.receivers.info")
611  do irec = 1,ig_nreceiver_quad
612  write(unit=100,fmt='(/,a,i10 )') "Receiver ",tg_receiver_quad(irec)%rglo
613  write(unit=100,fmt='( a,i10 )') " in core ",tg_receiver_quad(irec)%cpu-1
614  write(unit=100,fmt='( a,i10 )') " in quad ",tg_receiver_quad(irec)%iel
615  write(unit=100,fmt='( 2a )') " X "," Y "
616  write(unit=100,fmt='(2(e14.7,1x))') tg_receiver_quad(irec)%x,tg_receiver_quad(irec)%y
617  write(unit=100,fmt='( 2a )') " XI "," ETA "
618  write(unit=100,fmt='(2(e14.7,1x))') tg_receiver_quad(irec)%xi,tg_receiver_quad(irec)%eta
619  enddo
620  close(100)
621  endif
622 
623  else
624 
625  if (ig_myrank == 0) then
626  write(ig_lst_unit,'(/a)') "no free surface receiver computed"
627  call flush(ig_lst_unit)
628  endif
629 
630  endif
631 
632  return
633 !***********************************************************************************************************************************************************************************
634  end subroutine init_quad_receiver
635 !***********************************************************************************************************************************************************************************
636 
637 !
638 !
649 !***********************************************************************************************************************************************************************************
650  subroutine search_closest_hexa_gll(x,y,z,dmin,kgll,lgll,mgll,iel)
651 !***********************************************************************************************************************************************************************************
652 
653  use mpi
654 
655  use mod_global_variables, only :&
656  ig_ngll&
657  ,rg_gll_coordinate&
658  ,ig_nhexa&
659  ,ig_hexa_gll_glonum
660 
661  implicit none
662 
663  real , intent(in) :: x
664  real , intent(in) :: y
665  real , intent(in) :: z
666  real , intent(out) :: dmin
667  integer, intent(out) :: kgll
668  integer, intent(out) :: lgll
669  integer, intent(out) :: mgll
670  integer, intent(out) :: iel
671 
672  real :: d
673  integer :: ihexa
674  integer :: k
675  integer :: l
676  integer :: m
677  integer :: igll
678 
679  dmin = huge(dmin)
680 
681  do ihexa = 1,ig_nhexa
682  do k = 1,ig_ngll
683  do l = 1,ig_ngll
684  do m = 1,ig_ngll
685 
686  igll = ig_hexa_gll_glonum(m,l,k,ihexa)
687 
688  d = sqrt( (rg_gll_coordinate(1,igll) - x)**2 &
689  +(rg_gll_coordinate(2,igll) - y)**2 &
690  +(rg_gll_coordinate(3,igll) - z)**2 )
691 
692  if (d < dmin) then
693  dmin = d
694  kgll = k
695  lgll = l
696  mgll = m
697  iel = ihexa
698  endif
699 
700  enddo
701  enddo
702  enddo
703  enddo
704 
705  return
706 !***********************************************************************************************************************************************************************************
707  end subroutine search_closest_hexa_gll
708 !***********************************************************************************************************************************************************************************
709 
710 !
711 !
721 !***********************************************************************************************************************************************************************************
722  subroutine search_closest_quad_gll(x,y,global_gll,dmin,lgll,mgll,iel)
723 !***********************************************************************************************************************************************************************************
724 
725  use mpi
726 
727  use mod_global_variables, only :&
728  ig_ngll&
729  ,rg_gll_coordinate
730 
731  implicit none
732 
733  real , intent(in) :: x
734  real , intent(in) :: y
735  integer, intent(in), dimension(:,:,:) :: global_gll
736 
737  real , intent(out) :: dmin
738  integer, intent(out) :: lgll
739  integer, intent(out) :: mgll
740  integer, intent(out) :: iel
741 
742  real :: d
743 
744  integer :: iquad
745  integer :: l
746  integer :: m
747  integer :: nquad
748  integer :: igll
749 
750  dmin = huge(dmin)
751  nquad = size(global_gll,3)
752 
753  do iquad = 1,nquad
754  do l = 1,ig_ngll
755  do m = 1,ig_ngll
756 
757  igll = global_gll(m,l,iquad)
758  d = sqrt( (rg_gll_coordinate(1,igll) - x)**2 + (rg_gll_coordinate(2,igll) - y)**2 )
759 
760  if (d < dmin) then
761  dmin = d
762  lgll = l
763  mgll = m
764  iel = iquad
765  endif
766 
767  enddo
768  enddo
769  enddo
770 
771  return
772 !***********************************************************************************************************************************************************************************
773  end subroutine search_closest_quad_gll
774 !***********************************************************************************************************************************************************************************
775 
776 !
777 !
782 !***********************************************************************************************************************************************************************************
783  subroutine compute_info_hexa_receiver(tlxinf)
784 !***********************************************************************************************************************************************************************************
785 
786  use mpi
787 
788  use mod_global_variables, only :&
789  rg_gll_abscissa&
790  ,ig_ngll&
791  ,ig_hexa_gll_glonum&
792  ,type_receiver_hexa
793 
794  use mod_jacobian , only : compute_hexa_jacobian
795 
796  use mod_lagrange , only : lagrap
797 
798  use mod_coordinate , only : compute_hexa_point_coord
799 
800  implicit none
801 
802  real :: eps
803  real :: dmin
804  real :: xisol
805  real :: etsol
806  real :: zesol
807  real :: newx
808  real :: newy
809  real :: newz
810  real :: dxidx
811  real :: dxidy
812  real :: dxidz
813  real :: detdx
814  real :: detdy
815  real :: detdz
816  real :: dzedx
817  real :: dzedy
818  real :: dzedz
819  real :: dx
820  real :: dy
821  real :: dz
822  real :: dxi
823  real :: det
824  real :: dze
825  real :: deter
826 
827  integer :: ihexa
828  integer :: igll
829  integer :: k
830  integer :: l
831  integer :: m
832  integer :: iter
833 
834  integer, parameter :: iter_max=100
835 
836  type(type_receiver_hexa), intent(inout) :: tlxinf
837 
838 !
839 !---->initialize element number and coordinates (needed if dmin < eps)
840  ihexa= tlxinf%iel
841  newx = tlxinf%x
842  newy = tlxinf%y
843  newz = tlxinf%z
844  eps = 1.0e-3
845  iter = 0
846 !
847 !---->solve the nonlinear system to find local coordinates
848  xisol = rg_gll_abscissa(tlxinf%mgll)
849  etsol = rg_gll_abscissa(tlxinf%lgll)
850  zesol = rg_gll_abscissa(tlxinf%kgll)
851  dmin = tlxinf%dmin
852 
853  do while(dmin > eps)
854 
855  call compute_hexa_jacobian(ihexa,xisol,etsol,zesol,dxidx,dxidy,dxidz,detdx,detdy,detdz,dzedx,dzedy,dzedz,deter)
856 
857  call compute_hexa_point_coord(ihexa,xisol,etsol,zesol,newx,newy,newz)
858 
859  dx = tlxinf%x - newx
860  dy = tlxinf%y - newy
861  dz = tlxinf%z - newz
862  dxi = dxidx*dx + dxidy*dy + dxidz*dz
863  det = detdx*dx + detdy*dy + detdz*dz
864  dze = dzedx*dx + dzedy*dy + dzedz*dz
865  xisol = xisol + dxi
866  etsol = etsol + det
867  zesol = zesol + dze
868  dmin = sqrt( (newx-tlxinf%x)**2 + (newy-tlxinf%y)**2 + (newz-tlxinf%z)**2 )
869 
870  iter = iter + 1
871 
872  if (mod(iter,iter_max) == 0) then
873  eps = eps*2.0
874  endif
875 
876  enddo
877 
878  tlxinf%xi = xisol
879  tlxinf%eta = etsol
880  tlxinf%zeta = zesol
881  tlxinf%x = newx
882  tlxinf%y = newy
883  tlxinf%z = newz
884  tlxinf%dmin = dmin
885 !
886 !---->store global GLL nodes number
887  do k = 1,ig_ngll
888  do l = 1,ig_ngll
889  do m = 1,ig_ngll
890  igll = ig_hexa_gll_glonum(m,l,k,ihexa)
891  tlxinf%gll(m,l,k) = igll
892  enddo
893  enddo
894  enddo
895 !
896 !---->compute and store Lagrange polynomial values of a receiver at xisol, etsol and zesol
897  do k = 1,ig_ngll
898  do l = 1,ig_ngll
899  do m = 1,ig_ngll
900  tlxinf%lag(m,l,k) = lagrap(m,tlxinf%xi,ig_ngll)*lagrap(l,tlxinf%eta,ig_ngll)*lagrap(k,tlxinf%zeta,ig_ngll)
901  enddo
902  enddo
903  enddo
904 
905  return
906 !***********************************************************************************************************************************************************************************
907  end subroutine compute_info_hexa_receiver
908 !***********************************************************************************************************************************************************************************
909 
910 !
911 !
916 !***********************************************************************************************************************************************************************************
917  subroutine compute_info_quad_receiver(tlxinf)
918 !***********************************************************************************************************************************************************************************
919 
920  use mpi
921 
922  use mod_global_variables, only : rg_gll_abscissa&
923  ,ig_ngll&
924  ,ig_quadf_gll_glonum&
925  ,ig_myrank&
926  ,type_receiver_quad
927 
928  use mod_jacobian , only : compute_quad_jacobian
929 
930  use mod_lagrange , only : lagrap
931 
932  use mod_coordinate , only :&
933  compute_quad_point_coord&
934  ,compute_quad_point_coord_z
935 
936  implicit none
937 
938  type(type_receiver_quad), intent(inout) :: tlxinf
939 
940  real :: eps
941  real :: dmin
942  real :: xisol
943  real :: etsol
944  real :: newx
945  real :: newy
946  real :: newz
947  real :: dxidx
948  real :: dxidy
949  real :: detdx
950  real :: detdy
951  real :: deter
952  real :: dx
953  real :: dy
954  real :: dxi
955  real :: det
956 
957  integer :: k
958  integer :: l
959  integer :: iquad
960  integer :: igll
961  integer :: iter
962  integer, parameter :: iter_max=100
963 
964 !
965 !---->initialize coordinates, epsilon and iteration
966  newx = tlxinf%x
967  newy = tlxinf%y
968  dmin = tlxinf%dmin
969  iquad = tlxinf%iel
970  xisol = rg_gll_abscissa(tlxinf%mgll)
971  etsol = rg_gll_abscissa(tlxinf%lgll)
972  eps = 1.0e-2
973  iter = 0
974 !
975 !---->solve the nonlinear system to find local coordinates
976  do while(dmin > eps)
977 
978  call compute_quad_jacobian(iquad,xisol,etsol,dxidx,dxidy,detdx,detdy,deter)
979 
980  call compute_quad_point_coord(iquad,xisol,etsol,newx,newy,newz)
981 
982  dx = tlxinf%x - newx
983  dy = tlxinf%y - newy
984  dxi = dxidx*dx + dxidy*dy
985  det = detdx*dx + detdy*dy
986  xisol = xisol + dxi
987  etsol = etsol + det
988  dmin = sqrt( (newx-tlxinf%x)**2 + (newy-tlxinf%y)**2 )
989 
990  iter = iter + 1
991  if (mod(iter,iter_max) == 0) then
992  eps = eps*2.0
993  endif
994 
995  enddo
996 
997  tlxinf%xi = xisol
998  tlxinf%eta = etsol
999  tlxinf%x = newx
1000  tlxinf%y = newy
1001  tlxinf%z = compute_quad_point_coord_z(iquad,xisol,etsol)
1002  tlxinf%dmin = dmin
1003 !
1004 !---->store global GLL nodes number
1005  do k = 1,ig_ngll
1006  do l = 1,ig_ngll
1007  igll = ig_quadf_gll_glonum(l,k,iquad)
1008  tlxinf%gll(l,k) = igll
1009  enddo
1010  enddo
1011 !
1012 !---->compute and store Lagrange polynomial values of a receiver at xisol and etsol
1013  do k = 1,ig_ngll
1014  do l = 1,ig_ngll
1015  tlxinf%lag(l,k) = lagrap(l,tlxinf%xi,ig_ngll)*lagrap(k,tlxinf%eta,ig_ngll)
1016  enddo
1017  enddo
1018 
1019  return
1020 !***********************************************************************************************************************************************************************************
1021  end subroutine compute_info_quad_receiver
1022 !***********************************************************************************************************************************************************************************
1023 
1024 !
1025 !
1028 !***********************************************************************************************************************************************************************************
1029  subroutine write_receiver_output()
1030 !***********************************************************************************************************************************************************************************
1031 
1032  use mod_global_variables, only :&
1033  ig_nreceiver_hexa&
1034  ,ig_nreceiver_quad&
1035  ,rg_simu_current_time&
1036  ,tg_receiver_hexa&
1037  ,tg_receiver_quad&
1038  ,rg_gll_displacement&
1039  ,rg_gll_velocity&
1040  ,rg_gll_acceleration&
1041  ,ig_ngll&
1042  ,ig_ndof&
1043  ,ig_hexa_receiver_unit&
1044  ,ig_quad_receiver_unit
1045 
1046  use mod_gll_value
1047 
1048  use mod_lagrange, only :&
1049  hexa_lagrange_interp&
1050  ,quad_lagrange_interp
1051 
1052  implicit none
1053 
1054  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: gll_dis_hexa
1055  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: gll_vel_hexa
1056  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL,IG_NGLL) :: gll_acc_hexa
1057 
1058  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL) :: gll_dis_quad
1059  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL) :: gll_vel_quad
1060  real, dimension(IG_NDOF,IG_NGLL,IG_NGLL) :: gll_acc_quad
1061 
1062  real :: dis_x
1063  real :: dis_y
1064  real :: dis_z
1065  real :: vel_x
1066  real :: vel_y
1067  real :: vel_z
1068  real :: acc_x
1069  real :: acc_y
1070  real :: acc_z
1071 
1072  integer :: irec
1073 
1074 !
1075 !
1076 !********************************************************
1077 !---->receivers inside hexahedron elements
1078 !********************************************************
1079  do irec = 1,ig_nreceiver_hexa
1080 
1081  call get_hexa_gll_value(rg_gll_displacement,tg_receiver_hexa(irec)%gll,gll_dis_hexa)
1082  call get_hexa_gll_value(rg_gll_velocity ,tg_receiver_hexa(irec)%gll,gll_vel_hexa)
1083  call get_hexa_gll_value(rg_gll_acceleration,tg_receiver_hexa(irec)%gll,gll_acc_hexa)
1084 
1085 !
1086 !------->interpolate displacement, velocity and acceleration for receivers inside cpu myrank
1087  call hexa_lagrange_interp(gll_dis_hexa,tg_receiver_hexa(irec)%lag,dis_x,dis_y,dis_z)
1088  call hexa_lagrange_interp(gll_vel_hexa,tg_receiver_hexa(irec)%lag,vel_x,vel_y,vel_z)
1089  call hexa_lagrange_interp(gll_acc_hexa,tg_receiver_hexa(irec)%lag,acc_x,acc_y,acc_z)
1090 !
1091 !------->write binary values (files have been open inside the subroutine init_hexa_receiver)
1092  write(unit=ig_hexa_receiver_unit(irec)) rg_simu_current_time,dis_x,dis_y,dis_z,vel_x,vel_y,vel_z,acc_x,acc_y,acc_z
1093  call flush(ig_hexa_receiver_unit(irec))
1094 
1095  enddo
1096 
1097 !
1098 !
1099 !********************************************************
1100 !---->receivers inside quadrangle elements
1101 !********************************************************
1102  do irec = 1,ig_nreceiver_quad
1103 
1104  call get_quad_gll_value(rg_gll_displacement,tg_receiver_quad(irec)%gll,gll_dis_quad)
1105  call get_quad_gll_value(rg_gll_velocity ,tg_receiver_quad(irec)%gll,gll_vel_quad)
1106  call get_quad_gll_value(rg_gll_acceleration,tg_receiver_quad(irec)%gll,gll_acc_quad)
1107 
1108 !
1109 !------->interpolate displacement, velocity and acceleration for receivers inside cpu myrank
1110  call quad_lagrange_interp(gll_dis_quad,tg_receiver_quad(irec)%lag,dis_x,dis_y,dis_z)
1111  call quad_lagrange_interp(gll_vel_quad,tg_receiver_quad(irec)%lag,vel_x,vel_y,vel_z)
1112  call quad_lagrange_interp(gll_acc_quad,tg_receiver_quad(irec)%lag,acc_x,acc_y,acc_z)
1113 !
1114 !------->write binary values (files have been open inside the subroutine init_quad_receiver)
1115  write(unit=ig_quad_receiver_unit(irec)) rg_simu_current_time,dis_x,dis_y,dis_z,vel_x,vel_y,vel_z,acc_x,acc_y,acc_z
1116  call flush(ig_quad_receiver_unit(irec))
1117 
1118  enddo
1119 
1120  return
1121 
1122 !***********************************************************************************************************************************************************************************
1123  end subroutine write_receiver_output
1124 !***********************************************************************************************************************************************************************************
1125 
1126 
1127 !
1128 !
1130 !
1142 !***********************************************************************************************************************************************************************************
1143  subroutine is_inside_quad(x,y,quad_gnode,nquad,quad_nnode,dmin,lgll,mgll,quad_num,is_inside)
1144 !***********************************************************************************************************************************************************************************
1145 
1146  use mod_global_variables, only :&
1147  rg_gnode_x&
1148  ,rg_gnode_y
1149  implicit none
1150 
1151  real , intent( in) :: x
1152  real , intent( in) :: y
1153  integer, intent( in) :: nquad
1154  integer, intent( in) :: quad_nnode
1155  integer, intent( in), dimension(quad_nnode,nquad) :: quad_gnode
1156  real , intent(out) :: dmin
1157  integer, intent(out) :: lgll
1158  integer, intent(out) :: mgll
1159  integer, intent(out) :: quad_num
1160  logical, intent(out) :: is_inside
1161 
1162  integer :: node_num
1163  integer :: inode
1164  integer :: iquad
1165  integer :: l
1166  integer :: m
1167  real , dimension(quad_nnode+1) :: quad_gnode_x
1168  real , dimension(quad_nnode+1) :: quad_gnode_y
1169 
1170 
1171  is_inside = .false.
1172  dmin = huge(dmin)
1173 
1174  do iquad = 1,nquad
1175 
1176  do inode = 1,quad_nnode
1177 
1178  node_num = quad_gnode(inode,iquad)
1179 
1180  quad_gnode_x(inode) = rg_gnode_x(node_num)
1181  quad_gnode_y(inode) = rg_gnode_y(node_num)
1182 
1183  enddo
1184 
1185  quad_gnode_x(quad_nnode+1) = quad_gnode_x(1)
1186  quad_gnode_y(quad_nnode+1) = quad_gnode_y(1)
1187 
1188  call is_inside_polygon(x,y,quad_gnode_x,quad_gnode_y,quad_nnode+1,l,m)
1189 
1190  if (l >= 0) then
1191 
1192  is_inside = .true.
1193  quad_num = iquad
1194  lgll = 1
1195  mgll = 1
1196  dmin = sqrt( (quad_gnode_x(1)-x)**2 + (quad_gnode_y(1)-y)**2 )
1197 
1198  return
1199 
1200  endif
1201 
1202  enddo
1203 
1204  return
1205 
1206 !***********************************************************************************************************************************************************************************
1207  end subroutine is_inside_quad
1208 !**********************************************************************************************************************************************************************************
1209 
1210 
1211 !
1212 !
1214 !
1216 !
1218 !
1228 !***********************************************************************************************************************************************************************************
1229  subroutine is_inside_polygon(x0,y0,x,y,n,l,m)
1230 !***********************************************************************************************************************************************************************************
1231 ! Given a polygonal line connecting the vertices (x(i),y(i)) (i = 1,...,n) taken in this order.
1232 ! It is assumed that the polygonal path is a loop, where (x(n),y(n)) = (x(1),y(1)) or there is
1233 ! an arc from (x(n),y(n)) to (x(1),y(1)).
1234 !
1235 ! N.B.: the polygon may cross itself any number of times.
1236 !
1237 ! (x0,y0) is an arbitrary point and l and m are variables.
1238 ! on output, l and m are assigned the following values:
1239 ! l = -1 if (x0,y0) is outside the polygonal path
1240 ! l = 0 if (x0,y0) lies on the polygonal path
1241 ! l = 1 if (x0,y0) is inside the polygonal path
1242 ! m = 0 if (x0,y0) is on or outside the path. if (x0,y0) is inside the
1243 ! path then m is the winding number of the path around the point (x0,y0).
1244 !
1245  use mod_global_variables, only : rg_pi
1246 
1247  implicit none
1248 
1249  real , intent(in) :: x0
1250  real , intent(in) :: y0
1251  real , intent(in) :: x(n)
1252  real , intent(in) :: y(n)
1253  integer , intent(in) :: n
1254  integer , intent(out) :: l
1255  integer , intent(out) :: m
1256 
1257  real , parameter :: pix2 = 2.0*rg_pi
1258 
1259  integer :: i
1260  integer :: n0
1261  real :: angle
1262  real :: eps
1263  real :: sum
1264  real :: theta
1265  real :: theta1
1266  real :: thetai
1267  real :: tol
1268  real :: u
1269  real :: v
1270 
1271 
1272  eps = epsilon(eps)
1273 
1274  n0 = n
1275  if ( (x(1) == x(n)) .and. (y(1) == y(n)) ) n0 = n - 1
1276 
1277  tol = 4.0*eps*rg_pi
1278  l = -1
1279  m = 0
1280 
1281  u = x(1) - x0
1282  v = y(1) - y0
1283  if (u == 0.0 .and. v == 0.0) then
1284  l = 0
1285  return
1286  endif
1287 
1288  if (n0 < 2) return
1289  theta1 = atan2(v, u)
1290 
1291  sum = 0.0
1292  theta = theta1
1293 
1294  do i = 2,n0
1295 
1296  u = x(i) - x0
1297  v = y(i) - y0
1298  if (u == 0.0 .and. v == 0.0) then
1299  l = 0
1300  return
1301  endif
1302 
1303  thetai = atan2(v, u)
1304  angle = abs(thetai - theta)
1305 
1306  if (abs(angle - rg_pi) < tol) then
1307  l = 0
1308  return
1309  endif
1310 
1311  if (angle > rg_pi) angle = angle - pix2
1312  if (theta > thetai ) angle = -angle
1313 
1314  sum = sum + angle
1315  theta = thetai
1316 
1317  enddo
1318 
1319  angle = abs(theta1 - theta)
1320  if (abs(angle - rg_pi) < tol) then
1321  l = 0
1322  return
1323  endif
1324  if (angle > rg_pi ) angle = angle - pix2
1325  if (theta > theta1 ) angle = -angle
1326  sum = sum + angle
1327  !
1328  !sum = 2*RG_PI*m where m is the winding number
1329  m = abs(sum)/pix2 + 0.2
1330  if (m == 0) return
1331  l = 1
1332  if (sum < 0.0) m = -m
1333  return
1334 
1335 !***********************************************************************************************************************************************************************************
1336  end subroutine is_inside_polygon
1337 !***********************************************************************************************************************************************************************************
1338 
1339 end module mod_receiver
mod_receiver::search_closest_quad_gll
subroutine, public search_closest_quad_gll(x, y, global_gll, dmin, lgll, mgll, iel)
This subroutine searches among all free surface quadrangle elements in cpu myrank the closest GLL nod...
Definition: module_receiver.f90:722
mod_receiver
This module contains subroutines to compute information about receivers and to write receivers&#39; time ...
Definition: module_receiver.f90:126
mod_lagrange::lagrap
real function, public lagrap(i, x, n)
function to compute value of order Lagrange polynomial of the GLL node i at abscissa : ...
Definition: module_lagrange.f90:150
mod_receiver::search_closest_hexa_gll
subroutine, public search_closest_hexa_gll(x, y, z, dmin, kgll, lgll, mgll, iel)
This subroutine searches among all hexahedron elements in cpu myrank the closest GLL node to a point ...
Definition: module_receiver.f90:650
mod_jacobian::compute_hexa_jacobian
subroutine, public compute_hexa_jacobian(ihexa, xisol, etsol, zesol, dxidx, dxidy, dxidz, detdx, detdy, detdz, dzedx, dzedy, dzedz, det)
This subroutine computes jacobian matrix and its determinant at location in hexahedron element ihexa...
Definition: module_jacobian.f90:158
mod_receiver::is_inside_polygon
subroutine, private is_inside_polygon(x0, y0, x, y, n, l, m)
This subroutine tests if a point P with coordinates is inside a polygonal line.
Definition: module_receiver.f90:1229
mod_global_variables::type_receiver_quad
type for receivers (i.e., stations) located inside quadrangle elements
Definition: module_global_variables.f90:753
mod_gll_value::get_hexa_gll_value
subroutine, public get_hexa_gll_value(gll_all_values, gll_global_number, gll_selected_values)
subroutine to get hexahedron element GLL values from an array ordered by global GLL nodes numbering...
Definition: module_gll_value.f90:154
mod_global_variables
This module defines all global variables of EFISPEC3D. Scalar variables are initialized directly in t...
Definition: module_global_variables.f90:128
mod_receiver::init_quad_receiver
subroutine, public init_quad_receiver()
This subroutine reads all receivers in file *.fsr; determines to which cpu they belong; computes thei...
Definition: module_receiver.f90:371
mod_receiver::compute_info_quad_receiver
subroutine, public compute_info_quad_receiver(tlxinf)
This subroutine computes local coordinates of a receiver in the quadrangle element it belongs...
Definition: module_receiver.f90:917
mod_coordinate
This module contains subroutines to compute global -coordinates of a given point in the physical doma...
Definition: module_coordinate.f90:126
mod_coordinate::compute_hexa_point_coord
subroutine, public compute_hexa_point_coord(ihexa, xi, eta, zeta, x, y, z)
subroutine to compute -coordinates of a point knowing to which hexahedron element it belongs and its ...
Definition: module_coordinate.f90:151
mod_coordinate::compute_quad_point_coord_z
real function, public compute_quad_point_coord_z(iquad, xi, eta)
function to compute -coordinates of a point knowing to which quadrangle element it belongs and its lo...
Definition: module_coordinate.f90:287
mod_lagrange::quad_lagrange_interp
subroutine, public quad_lagrange_interp(gll_values, lag, x, y, z)
This subroutine interpolates GLL nodes -values of a quadrangle element using Lagrange polynomials...
Definition: module_lagrange.f90:365
mod_global_variables::type_receiver_hexa
type for receivers (i.e., stations) located inside hexahedron elements
Definition: module_global_variables.f90:726
mod_receiver::init_hexa_receiver
subroutine, public init_hexa_receiver()
This subroutine reads all receivers in file *.vor; determines to which cpu they belong; computes thei...
Definition: module_receiver.f90:154
mod_receiver::write_receiver_output
subroutine, public write_receiver_output()
This subroutine writes -displacements, velocities and accelerations at receivers. ...
Definition: module_receiver.f90:1029
mod_gll_value
This module contains subroutines to get GLL nodes values from global GLL nodes numbering.
Definition: module_gll_value.f90:126
mod_coordinate::compute_quad_point_coord
subroutine, public compute_quad_point_coord(iquad, xi, eta, x, y, z)
subroutine to compute -coordinates of a point knowing to which quadrangle element it belongs and its ...
Definition: module_coordinate.f90:223
mod_jacobian::compute_quad_jacobian
subroutine, public compute_quad_jacobian(iquad, xisol, etsol, dxidx, dxidy, detdx, detdy, det)
This subroutine computes jacobian matrix and its determinant at location in quadrangle element iquad...
Definition: module_jacobian.f90:305
mod_gll_value::get_quad_gll_value
subroutine, public get_quad_gll_value(gll_all_values, gll_global_number, gll_selected_values)
subroutine to get quadrangle element GLL values from an array ordered by global GLL nodes numbering...
Definition: module_gll_value.f90:200
mod_receiver::is_inside_quad
subroutine, public is_inside_quad(x, y, quad_gnode, nquad, quad_nnode, dmin, lgll, mgll, quad_num, is_inside)
This subroutine finds which quadrangle element contains a receiver.
Definition: module_receiver.f90:1143
mod_lagrange
This module contains functions to compute Lagrange polynomials and its derivatives; and to interpolat...
Definition: module_lagrange.f90:126
mod_receiver::compute_info_hexa_receiver
subroutine, public compute_info_hexa_receiver(tlxinf)
This subroutine computes local coordinates of a receiver in the hexahedron element it belongs...
Definition: module_receiver.f90:783
mod_jacobian
This module contains subroutines to compute jacobian matrix.
Definition: module_jacobian.f90:126
mod_lagrange::hexa_lagrange_interp
subroutine, public hexa_lagrange_interp(gll_values, lag, x, y, z)
This subroutine interpolates GLL nodes -values of a hexahedron element using Lagrange polynomials...
Definition: module_lagrange.f90:314