All Classes Files Functions Variables Pages
module_jacobian.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_jacobian
127 
128  use mpi
129 
130  implicit none
131 
132  private
133 
134  public :: compute_hexa_jacobian
135  public :: compute_quad_jacobian
136 
137  contains
138 
139 !
140 !
157 !***********************************************************************************************************************************************************************************
158  subroutine compute_hexa_jacobian(ihexa,xisol,etsol,zesol,dxidx,dxidy,dxidz,detdx,detdy,detdz,dzedx,dzedy,dzedz,det)
159 !***********************************************************************************************************************************************************************************
160 
161  use mpi
162 
163  use mod_global_variables, only : ig_hexa_gnode_xiloc&
164  ,ig_hexa_gnode_etloc&
165  ,ig_hexa_gnode_zeloc&
166  ,ig_line_nnode&
167  ,rg_gnode_x&
168  ,rg_gnode_y&
169  ,rg_gnode_z&
170  ,ig_hexa_nnode&
171  ,ig_hexa_gnode_glonum&
172  ,ig_myrank
173 
174  use mod_lagrange , only :&
175  lagrap_geom&
176  ,lagrad_geom
177 
178  implicit none
179 
180  real , intent( in) :: xisol
181  real , intent( in) :: etsol
182  real , intent( in) :: zesol
183  integer, intent( in) :: ihexa
184  real , intent(out) :: dxidx
185  real , intent(out) :: dxidy
186  real , intent(out) :: dxidz
187  real , intent(out) :: detdx
188  real , intent(out) :: detdy
189  real , intent(out) :: detdz
190  real , intent(out) :: dzedx
191  real , intent(out) :: dzedy
192  real , intent(out) :: dzedz
193  real , intent(out) :: det
194 
195  real , parameter :: eps = 1.0e-8
196  real :: inv_det
197  real :: xxi
198  real :: xet
199  real :: xze
200  real :: yxi
201  real :: yet
202  real :: yze
203  real :: zxi
204  real :: zet
205  real :: zze
206  real :: lag_deriv_xi
207  real :: lag_deriv_et
208  real :: lag_deriv_ze
209  real :: lag_xi
210  real :: lag_et
211  real :: lag_ze
212  real :: xgnode
213  real :: ygnode
214  real :: zgnode
215  integer :: m
216  integer :: ios
217 
218 !
219 !
220 !---->initialize partial derivatives
221  xxi = 0.0
222  xet = 0.0
223  xze = 0.0
224  yxi = 0.0
225  yet = 0.0
226  yze = 0.0
227  zxi = 0.0
228  zet = 0.0
229  zze = 0.0
230 
231 !
232 !
233 !---->compute partial derivatives
234  do m = 1,ig_hexa_nnode
235 
236  lag_deriv_xi = lagrad_geom(ig_hexa_gnode_xiloc(m),xisol,ig_line_nnode)
237  lag_deriv_et = lagrad_geom(ig_hexa_gnode_etloc(m),etsol,ig_line_nnode)
238  lag_deriv_ze = lagrad_geom(ig_hexa_gnode_zeloc(m),zesol,ig_line_nnode)
239 
240  lag_xi = lagrap_geom(ig_hexa_gnode_xiloc(m),xisol,ig_line_nnode)
241  lag_et = lagrap_geom(ig_hexa_gnode_etloc(m),etsol,ig_line_nnode)
242  lag_ze = lagrap_geom(ig_hexa_gnode_zeloc(m),zesol,ig_line_nnode)
243 
244  xgnode = rg_gnode_x(ig_hexa_gnode_glonum(m,ihexa))
245  ygnode = rg_gnode_y(ig_hexa_gnode_glonum(m,ihexa))
246  zgnode = rg_gnode_z(ig_hexa_gnode_glonum(m,ihexa))
247 
248  xxi = xxi + lag_deriv_xi*lag_et*lag_ze*xgnode
249  yxi = yxi + lag_deriv_xi*lag_et*lag_ze*ygnode
250  zxi = zxi + lag_deriv_xi*lag_et*lag_ze*zgnode
251 
252  xet = xet + lag_xi*lag_deriv_et*lag_ze*xgnode
253  yet = yet + lag_xi*lag_deriv_et*lag_ze*ygnode
254  zet = zet + lag_xi*lag_deriv_et*lag_ze*zgnode
255 
256  xze = xze + lag_xi*lag_et*lag_deriv_ze*xgnode
257  yze = yze + lag_xi*lag_et*lag_deriv_ze*ygnode
258  zze = zze + lag_xi*lag_et*lag_deriv_ze*zgnode
259 
260  enddo
261 !
262 !
263 !---->compute determinant and its inverse
264  det = xxi*yet*zze - xxi*yze*zet + xet*yze*zxi - xet*yxi*zze + xze*yxi*zet - xze*yet*zxi
265  inv_det = 1.0/det
266 
267  if (det.lt.eps) then
268  write(*,*) "***error in compute_hexa_jacobian: det null or negative, det = ",det
269  write(*,*) "***element ",ihexa, "cpu ",ig_myrank
270  call mpi_abort(mpi_comm_world,100,ios)
271  stop
272  endif
273 
274 !
275 !
276 !---->compute partial derivatives
277  dxidx = (yet*zze-yze*zet)*inv_det
278  dxidy = (xze*zet-xet*zze)*inv_det
279  dxidz = (xet*yze-xze*yet)*inv_det
280  detdx = (yze*zxi-yxi*zze)*inv_det
281  detdy = (xxi*zze-xze*zxi)*inv_det
282  detdz = (xze*yxi-xxi*yze)*inv_det
283  dzedx = (yxi*zet-yet*zxi)*inv_det
284  dzedy = (xet*zxi-xxi*zet)*inv_det
285  dzedz = (xxi*yet-xet*yxi)*inv_det
286 
287  return
288 !***********************************************************************************************************************************************************************************
289  end subroutine compute_hexa_jacobian
290 !***********************************************************************************************************************************************************************************
291 
292 !
293 !
304 !***********************************************************************************************************************************************************************************
305  subroutine compute_quad_jacobian(iquad,xisol,etsol,dxidx,dxidy,detdx,detdy,det)
306 !***********************************************************************************************************************************************************************************
307 
308  use mpi
309 
310  use mod_global_variables, only : ig_quad_gnode_xiloc&
311  ,ig_quad_gnode_etloc&
312  ,ig_line_nnode&
313  ,rg_gnode_x&
314  ,rg_gnode_y&
315  ,rg_gnode_z&
316  ,ig_quad_nnode&
317  ,ig_quadf_gnode_glonum&
318  ,ig_myrank
319 
320  use mod_lagrange , only :&
321  lagrap_geom&
322  ,lagrad_geom
323  implicit none
324 
325  real , intent(out) :: dxidx
326  real , intent(out) :: dxidy
327  real , intent(out) :: detdx
328  real , intent(out) :: detdy
329  real , intent( in) :: xisol
330  real , intent( in) :: etsol
331  integer, intent( in) :: iquad
332 
333  real , parameter :: eps = 1.0e-8
334  real :: xxi
335  real :: xet
336  real :: yxi
337  real :: yet
338  real :: det
339  real :: inv_det
340  real :: lag_deriv_xi
341  real :: lag_deriv_et
342  real :: lag_xi
343  real :: lag_et
344  real :: xgnode
345  real :: ygnode
346  integer :: m
347  integer :: ios
348 
349 !
350 !
351 !---->initialize partial derivatives
352  xxi = 0.0
353  xet = 0.0
354  yxi = 0.0
355  yet = 0.0
356 
357 !
358 !
359 !---->compute partial derivatives
360  do m = 1,ig_quad_nnode
361 
362  lag_deriv_xi = lagrad_geom(ig_quad_gnode_xiloc(m),xisol,ig_line_nnode)
363  lag_deriv_et = lagrad_geom(ig_quad_gnode_etloc(m),etsol,ig_line_nnode)
364 
365  lag_xi = lagrap_geom(ig_quad_gnode_xiloc(m),xisol,ig_line_nnode)
366  lag_et = lagrap_geom(ig_quad_gnode_etloc(m),etsol,ig_line_nnode)
367 
368  xgnode = rg_gnode_x(ig_quadf_gnode_glonum(m,iquad))
369  ygnode = rg_gnode_y(ig_quadf_gnode_glonum(m,iquad))
370 
371  xxi = xxi + lag_deriv_xi*lag_et*xgnode
372  yxi = yxi + lag_deriv_xi*lag_et*ygnode
373 
374  xet = xet + lag_xi*lag_deriv_et*xgnode
375  yet = yet + lag_xi*lag_deriv_et*ygnode
376 
377  enddo
378 
379 !
380 !
381 !---->compute determinant and its inverse
382  det = xxi*yet - xet*yxi
383  inv_det = 1.0/det
384 
385  if (det.lt.eps) then
386  write(*,*) "***error in compute_quad_jacobian: det null or negative, det = ",det
387  write(*,*) "***element ",iquad, "cpu ",ig_myrank
388  call mpi_abort(mpi_comm_world,100,ios)
389  stop
390  endif
391 
392 !
393 !
394 !---->compute partial derivatives
395  dxidx = (+yet)*inv_det
396  dxidy = (-xet)*inv_det
397  detdx = (-yxi)*inv_det
398  detdy = (+xxi)*inv_det
399 
400  return
401 !***********************************************************************************************************************************************************************************
402  end subroutine compute_quad_jacobian
403 !***********************************************************************************************************************************************************************************
404 
405 end module mod_jacobian
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_lagrange::lagrad_geom
real function, public lagrad_geom(i, x, n)
function to compute value of the derivative of order Lagrange polynomial of geometric node i at absc...
Definition: module_lagrange.f90:262
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_lagrange::lagrap_geom
real function, public lagrap_geom(i, x, n)
function to compute value of order Lagrange polynomial of geometric node i at abscissa : ...
Definition: module_lagrange.f90:181
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_lagrange
This module contains functions to compute Lagrange polynomials and its derivatives; and to interpolat...
Definition: module_lagrange.f90:126
mod_jacobian
This module contains subroutines to compute jacobian matrix.
Definition: module_jacobian.f90:126