diff --git a/utils/Visualization/opendx_AVS/create_OpenDX_files_to_view_a_3D_mesh.f90 b/utils/Visualization/opendx_AVS/create_OpenDX_files_to_view_a_3D_mesh.f90 new file mode 100644 index 000000000..cf0a70dd4 --- /dev/null +++ b/utils/Visualization/opendx_AVS/create_OpenDX_files_to_view_a_3D_mesh.f90 @@ -0,0 +1,173 @@ + +!===================================================================== +! +! S p e c f e m 3 D V e r s i o n 3 . 0 +! --------------------------------------- +! +! Main historical authors: Dimitri Komatitsch and Jeroen Tromp +! CNRS, France +! and Princeton University, USA +! (there are currently many more authors!) +! (c) October 2017 +! +! This program is free software; you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation; either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License along +! with this program; if not, write to the Free Software Foundation, Inc., +! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +! +!===================================================================== + + program create_OpenDX_file_to_view_the_mesh + +! Dimitri Komatitsch, CNRS Marseille, France, September 2018 + + implicit none + + integer :: nspec,nglob + integer :: ispec,iglob + integer :: iglob_read,ispec_loop,iformat,NGNOD,ia + + double precision, dimension(:), allocatable :: x,y,z + + double precision :: xread,yread,zread,xmin,xmax,ymin,ymax,zmin,zmax + + integer, dimension(:,:), allocatable :: ibool + integer, dimension(:), allocatable :: imaterial + +! read the mesh files in ASCII format (that is the standard case) + iformat = 1 + +! the mesh contains HEX8 elements + NGNOD = 8 + +! open SPECFEM3D_Cartesian mesh file to read the points + if (iformat == 1) then + open(unit=23,file='DATA_3D/nodes_coords_file',status='old',action='read') + read(23,*) nglob + else + open(unit=23,file='DATA_3D/nodes_coords_file.bin',form='unformatted',status='old',action='read') + read(23) nglob + endif + allocate(x(nglob)) + allocate(y(nglob)) + allocate(z(nglob)) + if (iformat == 1) then + do iglob = 1,nglob + read(23,*) iglob_read,xread,yread,zread + x(iglob_read) = xread + y(iglob_read) = yread + z(iglob_read) = zread + enddo + else + read(23) x + read(23) y + read(23) z + endif + close(23) + +! compute the min and max values of each coordinate + xmin = minval(x) + xmax = maxval(x) + + ymin = minval(y) + ymax = maxval(y) + + zmin = minval(z) + zmax = maxval(z) + + print *,'Xmin and Xmax of the mesh read = ',xmin,xmax + print *,'Ymin and Ymax of the mesh read = ',ymin,ymax + print *,'Zmin and Zmax of the mesh read = ',zmin,zmax + print * + +! ************* read mesh elements ************* + +! open SPECFEM3D_Cartesian topology file to read the mesh elements + if (iformat == 1) then + open(unit=23,file='DATA_3D/mesh_file',status='old',action='read') + read(23,*) nspec + else + open(unit=23,file='DATA_3D/mesh_file.bin',form='unformatted',status='old',action='read') + read(23) nspec + endif + + allocate(ibool(NGNOD,nspec)) + allocate(imaterial(nspec)) + +! loop on the whole mesh + if (iformat == 1) then + do ispec_loop = 1,nspec + read(23,*) ispec,(ibool(ia,ispec), ia = 1,NGNOD) + enddo + else + read(23) ibool + endif + + close(23) + + print *,'Total number of elements in the mesh read = ',nspec + print * + +! read the material property numbers + open(unit=23,file='DATA_3D/materials_file',status='old',action='read') + do ispec_loop = 1,nspec + read(23,*) ispec,imaterial(ispec) + enddo + close(23) + +!! DK DK create an OpenDX file showing all the elements of the mesh + + open(unit=11,file='DX_fullmesh.dx',status='unknown') + write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows' + + do iglob = 1,nglob + write(11,"(f10.7,1x,f10.7,1x,f10.7)") x(iglob),y(iglob),z(iglob) + enddo + +! ************* generate elements ****************** + + write(11,*) 'object 2 class array type int rank 1 shape 8 items ',nspec,' data follows' + +! read local elements in this slice and output global AVS or DX elements + do ispec=1,nspec +! point order in OpenDX in 2D is 1,4,2,3 *not* 1,2,3,4 as in AVS +! point order in OpenDX in 3D is 4,1,8,5,3,2,7,6, *not* 1,2,3,4,5,6,7,8 as in AVS +! in the case of OpenDX, node numbers start at zero + write(11,"(i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9)") & + ibool(4,ispec)-1, ibool(1,ispec)-1, ibool(8,ispec)-1, ibool(5,ispec)-1, & + ibool(3,ispec)-1, ibool(2,ispec)-1, ibool(7,ispec)-1, ibool(6,ispec)-1 + enddo + +! ************* generate element data values ****************** + +! output AVS or DX header for data + write(11,*) 'attribute "element type" string "cubes"' + write(11,*) 'attribute "ref" string "positions"' + write(11,*) 'object 3 class array type float rank 0 items ',nspec,' data follows' + +! read local elements in this slice and output global AVS or DX elements + do ispec=1,nspec + write(11,*) imaterial(ispec) + enddo + +! define OpenDX field + write(11,*) 'attribute "dep" string "connections"' + write(11,*) 'object "irregular positions irregular connections" class field' + write(11,*) 'component "positions" value 1' + write(11,*) 'component "connections" value 2' + write(11,*) 'component "data" value 3' + write(11,*) 'end' + + close(11) + + end program create_OpenDX_file_to_view_the_mesh + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_OpenDX_files_to_view_the_3D_mesh.f90 b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_OpenDX_files_to_view_the_3D_mesh.f90 new file mode 100644 index 000000000..cf0a70dd4 --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_OpenDX_files_to_view_the_3D_mesh.f90 @@ -0,0 +1,173 @@ + +!===================================================================== +! +! S p e c f e m 3 D V e r s i o n 3 . 0 +! --------------------------------------- +! +! Main historical authors: Dimitri Komatitsch and Jeroen Tromp +! CNRS, France +! and Princeton University, USA +! (there are currently many more authors!) +! (c) October 2017 +! +! This program is free software; you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation; either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License along +! with this program; if not, write to the Free Software Foundation, Inc., +! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +! +!===================================================================== + + program create_OpenDX_file_to_view_the_mesh + +! Dimitri Komatitsch, CNRS Marseille, France, September 2018 + + implicit none + + integer :: nspec,nglob + integer :: ispec,iglob + integer :: iglob_read,ispec_loop,iformat,NGNOD,ia + + double precision, dimension(:), allocatable :: x,y,z + + double precision :: xread,yread,zread,xmin,xmax,ymin,ymax,zmin,zmax + + integer, dimension(:,:), allocatable :: ibool + integer, dimension(:), allocatable :: imaterial + +! read the mesh files in ASCII format (that is the standard case) + iformat = 1 + +! the mesh contains HEX8 elements + NGNOD = 8 + +! open SPECFEM3D_Cartesian mesh file to read the points + if (iformat == 1) then + open(unit=23,file='DATA_3D/nodes_coords_file',status='old',action='read') + read(23,*) nglob + else + open(unit=23,file='DATA_3D/nodes_coords_file.bin',form='unformatted',status='old',action='read') + read(23) nglob + endif + allocate(x(nglob)) + allocate(y(nglob)) + allocate(z(nglob)) + if (iformat == 1) then + do iglob = 1,nglob + read(23,*) iglob_read,xread,yread,zread + x(iglob_read) = xread + y(iglob_read) = yread + z(iglob_read) = zread + enddo + else + read(23) x + read(23) y + read(23) z + endif + close(23) + +! compute the min and max values of each coordinate + xmin = minval(x) + xmax = maxval(x) + + ymin = minval(y) + ymax = maxval(y) + + zmin = minval(z) + zmax = maxval(z) + + print *,'Xmin and Xmax of the mesh read = ',xmin,xmax + print *,'Ymin and Ymax of the mesh read = ',ymin,ymax + print *,'Zmin and Zmax of the mesh read = ',zmin,zmax + print * + +! ************* read mesh elements ************* + +! open SPECFEM3D_Cartesian topology file to read the mesh elements + if (iformat == 1) then + open(unit=23,file='DATA_3D/mesh_file',status='old',action='read') + read(23,*) nspec + else + open(unit=23,file='DATA_3D/mesh_file.bin',form='unformatted',status='old',action='read') + read(23) nspec + endif + + allocate(ibool(NGNOD,nspec)) + allocate(imaterial(nspec)) + +! loop on the whole mesh + if (iformat == 1) then + do ispec_loop = 1,nspec + read(23,*) ispec,(ibool(ia,ispec), ia = 1,NGNOD) + enddo + else + read(23) ibool + endif + + close(23) + + print *,'Total number of elements in the mesh read = ',nspec + print * + +! read the material property numbers + open(unit=23,file='DATA_3D/materials_file',status='old',action='read') + do ispec_loop = 1,nspec + read(23,*) ispec,imaterial(ispec) + enddo + close(23) + +!! DK DK create an OpenDX file showing all the elements of the mesh + + open(unit=11,file='DX_fullmesh.dx',status='unknown') + write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows' + + do iglob = 1,nglob + write(11,"(f10.7,1x,f10.7,1x,f10.7)") x(iglob),y(iglob),z(iglob) + enddo + +! ************* generate elements ****************** + + write(11,*) 'object 2 class array type int rank 1 shape 8 items ',nspec,' data follows' + +! read local elements in this slice and output global AVS or DX elements + do ispec=1,nspec +! point order in OpenDX in 2D is 1,4,2,3 *not* 1,2,3,4 as in AVS +! point order in OpenDX in 3D is 4,1,8,5,3,2,7,6, *not* 1,2,3,4,5,6,7,8 as in AVS +! in the case of OpenDX, node numbers start at zero + write(11,"(i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9)") & + ibool(4,ispec)-1, ibool(1,ispec)-1, ibool(8,ispec)-1, ibool(5,ispec)-1, & + ibool(3,ispec)-1, ibool(2,ispec)-1, ibool(7,ispec)-1, ibool(6,ispec)-1 + enddo + +! ************* generate element data values ****************** + +! output AVS or DX header for data + write(11,*) 'attribute "element type" string "cubes"' + write(11,*) 'attribute "ref" string "positions"' + write(11,*) 'object 3 class array type float rank 0 items ',nspec,' data follows' + +! read local elements in this slice and output global AVS or DX elements + do ispec=1,nspec + write(11,*) imaterial(ispec) + enddo + +! define OpenDX field + write(11,*) 'attribute "dep" string "connections"' + write(11,*) 'object "irregular positions irregular connections" class field' + write(11,*) 'component "positions" value 1' + write(11,*) 'component "connections" value 2' + write(11,*) 'component "data" value 3' + write(11,*) 'end' + + close(11) + + end program create_OpenDX_file_to_view_the_mesh + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_database_files_for_external_faces_of_the_3D_mesh.f90 b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_database_files_for_external_faces_of_the_3D_mesh.f90 new file mode 100644 index 000000000..f082b5fec --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_database_files_for_external_faces_of_the_3D_mesh.f90 @@ -0,0 +1,1212 @@ +!===================================================================== +! +! S p e c f e m 3 D V e r s i o n 3 . 0 +! --------------------------------------- +! +! Main historical authors: Dimitri Komatitsch and Jeroen Tromp +! CNRS, France +! and Princeton University, USA +! (there are currently many more authors!) +! (c) October 2017 +! +! This program is free software; you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation; either version 3 of the License, or +! (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License along +! with this program; if not, write to the Free Software Foundation, Inc., +! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +! +!===================================================================== + + program create_database_files_for_external_faces_of_the_3D_mesh + +! Dimitri Komatitsch, CNRS Marseille, France, April 2013 and February 2017 and September 2018 + + implicit none + + integer :: nspec,npoin + integer :: ispec,ipoin + integer :: ipoin_read,ispec_loop,iformat,NGNOD,ia + integer :: i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13,i14,i15,i16,i17,i18,i19,i20,i21,i22,i23,i24,i25,i26,i27 + integer :: count_faces_found + + double precision, dimension(:), allocatable :: x,y,z + + double precision :: xread,yread,zread,xmin,xmax,ymin,ymax,zmin,zmax,limit,size_of_model + + logical :: already_found_a_face + + integer, dimension(:,:), allocatable :: ibool + +! to make sure coordinate roundoff problems do not occur, use a tolerance of 0.5% + double precision, parameter :: SMALL_PERCENTAGE_TOLERANCE = 1.005d0 + + double precision, parameter :: SMALL_RELATIVE_VALUE = 0.5d-3 + +! read the mesh files in ASCII format (that is the standard case) + iformat = 1 + +! the mesh contains HEX8 elements + NGNOD = 8 + +! open SPECFEM3D_Cartesian mesh file to read the points + if (iformat == 1) then + open(unit=23,file='DATA_3D/nodes_coords_file',status='old',action='read') + read(23,*) npoin + else + open(unit=23,file='DATA_3D/nodes_coords_file.bin',form='unformatted',status='old',action='read') + read(23) npoin + endif + allocate(x(npoin)) + allocate(y(npoin)) + allocate(z(npoin)) + if (iformat == 1) then + do ipoin = 1,npoin + read(23,*) ipoin_read,xread,yread,zread + x(ipoin_read) = xread + y(ipoin_read) = yread + z(ipoin_read) = zread + enddo + else + read(23) x + read(23) y + read(23) z + endif + close(23) + +! compute the min and max values of each coordinate + xmin = minval(x) + xmax = maxval(x) + + ymin = minval(y) + ymax = maxval(y) + + zmin = minval(z) + zmax = maxval(z) + + print *,'Xmin and Xmax of the mesh read = ',xmin,xmax + print *,'Ymin and Ymax of the mesh read = ',ymin,ymax + print *,'Zmin and Zmax of the mesh read = ',zmin,zmax + print * + +! ************* read mesh elements ************* + +! open SPECFEM3D_Cartesian topology file to read the mesh elements + if (iformat == 1) then + open(unit=23,file='DATA_3D/mesh_file',status='old',action='read') + read(23,*) nspec + else + open(unit=23,file='DATA_3D/mesh_file.bin',form='unformatted',status='old',action='read') + read(23) nspec + endif + + allocate(ibool(NGNOD,nspec)) + +! loop on the whole mesh + if (iformat == 1) then + do ispec_loop = 1,nspec + read(23,*) ispec,(ibool(ia,ispec), ia = 1,NGNOD) + enddo + else + read(23) ibool + endif + + close(23) + + print *,'Total number of elements in the mesh read = ',nspec + print * + +! ************* generate "absorbing_surface_file_xmin" ************* + +! first count the number of faces that are along that edge + + count_faces_found = 0 + +! Xmin + size_of_model = xmax - xmin + limit = xmin + SMALL_RELATIVE_VALUE*size_of_model + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (.true.) then + + already_found_a_face = .false. + +! test face 1 (bottom) + if (x(i1) < limit .and. x(i2) < limit .and. x(i3) < limit .and. x(i4) < limit) then + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 2 (top) + if (x(i5) < limit .and. x(i6) < limit .and. x(i7) < limit .and. x(i8) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 3 (left) + if (x(i1) < limit .and. x(i4) < limit .and. x(i8) < limit .and. x(i5) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 4 (right) + if (x(i2) < limit .and. x(i3) < limit .and. x(i7) < limit .and. x(i6) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 5 (front) + if (x(i1) < limit .and. x(i2) < limit .and. x(i6) < limit .and. x(i5) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 6 (back) + if (x(i4) < limit .and. x(i3) < limit .and. x(i7) < limit .and. x(i8) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + + endif + + enddo + + print *,'found ',count_faces_found,' full faces on face Xmin' + +!----------------------------- + +! open SPECFEM3D_Cartesian topology file to read the mesh elements + open(unit=24,file='DATA_3D/absorbing_surface_file_xmin',status='unknown',action='write') + +! write the total number of face elements + write(24,*) count_faces_found + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (NGNOD == 27) then + i9 = ibool(9,ispec) + i10 = ibool(10,ispec) + i11 = ibool(11,ispec) + i12 = ibool(12,ispec) + i13 = ibool(13,ispec) + i14 = ibool(14,ispec) + i15 = ibool(15,ispec) + i16 = ibool(16,ispec) + i17 = ibool(17,ispec) + i18 = ibool(18,ispec) + i19 = ibool(19,ispec) + i20 = ibool(20,ispec) + i21 = ibool(21,ispec) + i22 = ibool(22,ispec) + i23 = ibool(23,ispec) + i24 = ibool(24,ispec) + i25 = ibool(25,ispec) + i26 = ibool(26,ispec) + i27 = ibool(27,ispec) + endif + + if (.true.) then + +! for the six faces below it is important to make sure we write the four points +! in an order for which the normal to the face points outwards + +! test face 1 (bottom) + if (x(i1) < limit .and. x(i2) < limit .and. x(i3) < limit .and. x(i4) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i4,i3,i2,i1 + else + write(24,"(10(i9,1x))") ispec,i4,i3,i2,i1,i11,i10,i9,i12,i21 + endif + endif + +! test face 2 (top) + if (x(i5) < limit .and. x(i6) < limit .and. x(i7) < limit .and. x(i8) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i5,i6,i7,i8 + else + write(24,"(10(i9,1x))") ispec,i5,i6,i7,i8,i17,i18,i19,i20,i26 + endif + endif + +! test face 3 (left) + if (x(i1) < limit .and. x(i4) < limit .and. x(i8) < limit .and. x(i5) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i5,i8,i4 + else + write(24,"(10(i9,1x))") ispec,i1,i5,i8,i4,i13,i20,i16,i12,i25 + endif + endif + +! test face 4 (right) + if (x(i2) < limit .and. x(i3) < limit .and. x(i7) < limit .and. x(i6) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i2,i3,i7,i6 + else + write(24,"(10(i9,1x))") ispec,i2,i3,i7,i6,i10,i15,i18,i14,i23 + endif + endif + +! test face 5 (front) + if (x(i1) < limit .and. x(i2) < limit .and. x(i6) < limit .and. x(i5) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i2,i6,i5 + else + write(24,"(10(i9,1x))") ispec,i1,i2,i6,i5,i9,i14,i17,i13,i22 + endif + endif + +! test face 6 (back) + if (x(i4) < limit .and. x(i3) < limit .and. x(i7) < limit .and. x(i8) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i3,i4,i8,i7 + else + write(24,"(10(i9,1x))") ispec,i3,i4,i8,i7,i11,i16,i19,i15,i24 + endif + endif + + endif + + enddo + + close(24) + + print *,'File "absorbing_surface_file_xmin" has been successfully created' + print * + +! ************* generate "absorbing_surface_file_xmax" ************* + +! first count the number of faces that are along that edge + + count_faces_found = 0 + +! Xmax + size_of_model = xmax - xmin + limit = xmax - SMALL_RELATIVE_VALUE*size_of_model + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (.true.) then + + already_found_a_face = .false. + +! test face 1 (bottom) + if (x(i1) > limit .and. x(i2) > limit .and. x(i3) > limit .and. x(i4) > limit) then + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 2 (top) + if (x(i5) > limit .and. x(i6) > limit .and. x(i7) > limit .and. x(i8) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 3 (left) + if (x(i1) > limit .and. x(i4) > limit .and. x(i8) > limit .and. x(i5) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 4 (right) + if (x(i2) > limit .and. x(i3) > limit .and. x(i7) > limit .and. x(i6) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 5 (front) + if (x(i1) > limit .and. x(i2) > limit .and. x(i6) > limit .and. x(i5) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 6 (back) + if (x(i4) > limit .and. x(i3) > limit .and. x(i7) > limit .and. x(i8) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + + endif + + enddo + + print *,'found ',count_faces_found,' full faces on face Xmax' + +!----------------------------- + + open(unit=24,file='DATA_3D/absorbing_surface_file_xmax',status='unknown',action='write') + +! write the total number of face elements + write(24,*) count_faces_found + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (NGNOD == 27) then + i9 = ibool(9,ispec) + i10 = ibool(10,ispec) + i11 = ibool(11,ispec) + i12 = ibool(12,ispec) + i13 = ibool(13,ispec) + i14 = ibool(14,ispec) + i15 = ibool(15,ispec) + i16 = ibool(16,ispec) + i17 = ibool(17,ispec) + i18 = ibool(18,ispec) + i19 = ibool(19,ispec) + i20 = ibool(20,ispec) + i21 = ibool(21,ispec) + i22 = ibool(22,ispec) + i23 = ibool(23,ispec) + i24 = ibool(24,ispec) + i25 = ibool(25,ispec) + i26 = ibool(26,ispec) + i27 = ibool(27,ispec) + endif + + if (.true.) then + +! for the six faces below it is important to make sure we write the four points +! in an order for which the normal to the face points outwards + +! test face 1 (bottom) + if (x(i1) > limit .and. x(i2) > limit .and. x(i3) > limit .and. x(i4) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i4,i3,i2,i1 + else + write(24,"(10(i9,1x))") ispec,i4,i3,i2,i1,i11,i10,i9,i12,i21 + endif + endif + +! test face 2 (top) + if (x(i5) > limit .and. x(i6) > limit .and. x(i7) > limit .and. x(i8) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i5,i6,i7,i8 + else + write(24,"(10(i9,1x))") ispec,i5,i6,i7,i8,i17,i18,i19,i20,i26 + endif + endif + +! test face 3 (left) + if (x(i1) > limit .and. x(i4) > limit .and. x(i8) > limit .and. x(i5) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i5,i8,i4 + else + write(24,"(10(i9,1x))") ispec,i1,i5,i8,i4,i13,i20,i16,i12,i25 + endif + endif + +! test face 4 (right) + if (x(i2) > limit .and. x(i3) > limit .and. x(i7) > limit .and. x(i6) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i2,i3,i7,i6 + else + write(24,"(10(i9,1x))") ispec,i2,i3,i7,i6,i10,i15,i18,i14,i23 + endif + endif + +! test face 5 (front) + if (x(i1) > limit .and. x(i2) > limit .and. x(i6) > limit .and. x(i5) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i2,i6,i5 + else + write(24,"(10(i9,1x))") ispec,i1,i2,i6,i5,i9,i14,i17,i13,i22 + endif + endif + +! test face 6 (back) + if (x(i4) > limit .and. x(i3) > limit .and. x(i7) > limit .and. x(i8) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i3,i4,i8,i7 + else + write(24,"(10(i9,1x))") ispec,i3,i4,i8,i7,i11,i16,i19,i15,i24 + endif + endif + + endif + + enddo + + close(24) + + print *,'File "absorbing_surface_file_xmax" has been successfully created' + print * + +! ************* generate "absorbing_surface_file_ymin" ************* + +! first count the number of faces that are along that edge + + count_faces_found = 0 + +! Ymin + size_of_model = ymax - ymin + limit = ymin + SMALL_RELATIVE_VALUE*size_of_model + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (.true.) then + + already_found_a_face = .false. + +! test face 1 (bottom) + if (y(i1) < limit .and. y(i2) < limit .and. y(i3) < limit .and. y(i4) < limit) then + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 2 (top) + if (y(i5) < limit .and. y(i6) < limit .and. y(i7) < limit .and. y(i8) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 3 (left) + if (y(i1) < limit .and. y(i4) < limit .and. y(i8) < limit .and. y(i5) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 4 (right) + if (y(i2) < limit .and. y(i3) < limit .and. y(i7) < limit .and. y(i6) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 5 (front) + if (y(i1) < limit .and. y(i2) < limit .and. y(i6) < limit .and. y(i5) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 6 (back) + if (y(i4) < limit .and. y(i3) < limit .and. y(i7) < limit .and. y(i8) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + + endif + + enddo + + print *,'found ',count_faces_found,' full faces on face Ymin' + +!----------------------------- + + open(unit=24,file='DATA_3D/absorbing_surface_file_ymin',status='unknown',action='write') + +! write the total number of face elements + write(24,*) count_faces_found + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (NGNOD == 27) then + i9 = ibool(9,ispec) + i10 = ibool(10,ispec) + i11 = ibool(11,ispec) + i12 = ibool(12,ispec) + i13 = ibool(13,ispec) + i14 = ibool(14,ispec) + i15 = ibool(15,ispec) + i16 = ibool(16,ispec) + i17 = ibool(17,ispec) + i18 = ibool(18,ispec) + i19 = ibool(19,ispec) + i20 = ibool(20,ispec) + i21 = ibool(21,ispec) + i22 = ibool(22,ispec) + i23 = ibool(23,ispec) + i24 = ibool(24,ispec) + i25 = ibool(25,ispec) + i26 = ibool(26,ispec) + i27 = ibool(27,ispec) + endif + + if (.true.) then + +! for the six faces below it is important to make sure we write the four points +! in an order for which the normal to the face points outwards + +! test face 1 (bottom) + if (y(i1) < limit .and. y(i2) < limit .and. y(i3) < limit .and. y(i4) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i4,i3,i2,i1 + else + write(24,"(10(i9,1x))") ispec,i4,i3,i2,i1,i11,i10,i9,i12,i21 + endif + endif + +! test face 2 (top) + if (y(i5) < limit .and. y(i6) < limit .and. y(i7) < limit .and. y(i8) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i5,i6,i7,i8 + else + write(24,"(10(i9,1x))") ispec,i5,i6,i7,i8,i17,i18,i19,i20,i26 + endif + endif + +! test face 3 (left) + if (y(i1) < limit .and. y(i4) < limit .and. y(i8) < limit .and. y(i5) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i5,i8,i4 + else + write(24,"(10(i9,1x))") ispec,i1,i5,i8,i4,i13,i20,i16,i12,i25 + endif + endif + +! test face 4 (right) + if (y(i2) < limit .and. y(i3) < limit .and. y(i7) < limit .and. y(i6) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i2,i3,i7,i6 + else + write(24,"(10(i9,1x))") ispec,i2,i3,i7,i6,i10,i15,i18,i14,i23 + endif + endif + +! test face 5 (front) + if (y(i1) < limit .and. y(i2) < limit .and. y(i6) < limit .and. y(i5) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i2,i6,i5 + else + write(24,"(10(i9,1x))") ispec,i1,i2,i6,i5,i9,i14,i17,i13,i22 + endif + endif + +! test face 6 (back) + if (y(i4) < limit .and. y(i3) < limit .and. y(i7) < limit .and. y(i8) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i3,i4,i8,i7 + else + write(24,"(10(i9,1x))") ispec,i3,i4,i8,i7,i11,i16,i19,i15,i24 + endif + endif + + endif + + enddo + + close(24) + + print *,'File "absorbing_surface_file_ymin" has been successfully created' + print * + +! ************* generate "absorbing_surface_file_ymax" ************* + +! first count the number of faces that are along that edge + + count_faces_found = 0 + +! Ymax + size_of_model = ymax - ymin + limit = ymax - SMALL_RELATIVE_VALUE*size_of_model + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (.true.) then + + already_found_a_face = .false. + +! test face 1 (bottom) + if (y(i1) > limit .and. y(i2) > limit .and. y(i3) > limit .and. y(i4) > limit) then + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 2 (top) + if (y(i5) > limit .and. y(i6) > limit .and. y(i7) > limit .and. y(i8) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 3 (left) + if (y(i1) > limit .and. y(i4) > limit .and. y(i8) > limit .and. y(i5) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 4 (right) + if (y(i2) > limit .and. y(i3) > limit .and. y(i7) > limit .and. y(i6) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 5 (front) + if (y(i1) > limit .and. y(i2) > limit .and. y(i6) > limit .and. y(i5) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 6 (back) + if (y(i4) > limit .and. y(i3) > limit .and. y(i7) > limit .and. y(i8) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + + endif + + enddo + + print *,'found ',count_faces_found,' full faces on face Ymax' + +!----------------------------- + + open(unit=24,file='DATA_3D/absorbing_surface_file_ymax',status='unknown',action='write') + +! write the total number of face elements + write(24,*) count_faces_found + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (NGNOD == 27) then + i9 = ibool(9,ispec) + i10 = ibool(10,ispec) + i11 = ibool(11,ispec) + i12 = ibool(12,ispec) + i13 = ibool(13,ispec) + i14 = ibool(14,ispec) + i15 = ibool(15,ispec) + i16 = ibool(16,ispec) + i17 = ibool(17,ispec) + i18 = ibool(18,ispec) + i19 = ibool(19,ispec) + i20 = ibool(20,ispec) + i21 = ibool(21,ispec) + i22 = ibool(22,ispec) + i23 = ibool(23,ispec) + i24 = ibool(24,ispec) + i25 = ibool(25,ispec) + i26 = ibool(26,ispec) + i27 = ibool(27,ispec) + endif + + if (.true.) then + +! for the six faces below it is important to make sure we write the four points +! in an order for which the normal to the face points outwards + +! test face 1 (bottom) + if (y(i1) > limit .and. y(i2) > limit .and. y(i3) > limit .and. y(i4) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i4,i3,i2,i1 + else + write(24,"(10(i9,1x))") ispec,i4,i3,i2,i1,i11,i10,i9,i12,i21 + endif + endif + +! test face 2 (top) + if (y(i5) > limit .and. y(i6) > limit .and. y(i7) > limit .and. y(i8) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i5,i6,i7,i8 + else + write(24,"(10(i9,1x))") ispec,i5,i6,i7,i8,i17,i18,i19,i20,i26 + endif + endif + +! test face 3 (left) + if (y(i1) > limit .and. y(i4) > limit .and. y(i8) > limit .and. y(i5) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i5,i8,i4 + else + write(24,"(10(i9,1x))") ispec,i1,i5,i8,i4,i13,i20,i16,i12,i25 + endif + endif + +! test face 4 (right) + if (y(i2) > limit .and. y(i3) > limit .and. y(i7) > limit .and. y(i6) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i2,i3,i7,i6 + else + write(24,"(10(i9,1x))") ispec,i2,i3,i7,i6,i10,i15,i18,i14,i23 + endif + endif + +! test face 5 (front) + if (y(i1) > limit .and. y(i2) > limit .and. y(i6) > limit .and. y(i5) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i2,i6,i5 + else + write(24,"(10(i9,1x))") ispec,i1,i2,i6,i5,i9,i14,i17,i13,i22 + endif + endif + +! test face 6 (back) + if (y(i4) > limit .and. y(i3) > limit .and. y(i7) > limit .and. y(i8) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i3,i4,i8,i7 + else + write(24,"(10(i9,1x))") ispec,i3,i4,i8,i7,i11,i16,i19,i15,i24 + endif + endif + + endif + + enddo + + close(24) + + print *,'File "absorbing_surface_file_ymax" has been successfully created' + print * + +! ************* generate "absorbing_surface_file_bottom" ************* + +! first count the number of faces that are along that edge + + count_faces_found = 0 + +! Zmin + size_of_model = zmax - zmin + limit = zmin + SMALL_RELATIVE_VALUE*size_of_model + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (.true.) then + + already_found_a_face = .false. + +! test face 1 (bottom) + if (z(i1) < limit .and. z(i2) < limit .and. z(i3) < limit .and. z(i4) < limit) then + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 2 (top) + if (z(i5) < limit .and. z(i6) < limit .and. z(i7) < limit .and. z(i8) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 3 (left) + if (z(i1) < limit .and. z(i4) < limit .and. z(i8) < limit .and. z(i5) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 4 (right) + if (z(i2) < limit .and. z(i3) < limit .and. z(i7) < limit .and. z(i6) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 5 (front) + if (z(i1) < limit .and. z(i2) < limit .and. z(i6) < limit .and. z(i5) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 6 (back) + if (z(i4) < limit .and. z(i3) < limit .and. z(i7) < limit .and. z(i8) < limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + + endif + + enddo + + print *,'found ',count_faces_found,' full faces on face Zmin' + +!----------------------------- + + open(unit=24,file='DATA_3D/absorbing_surface_file_bottom',status='unknown',action='write') + +! write the total number of face elements + write(24,*) count_faces_found + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (NGNOD == 27) then + i9 = ibool(9,ispec) + i10 = ibool(10,ispec) + i11 = ibool(11,ispec) + i12 = ibool(12,ispec) + i13 = ibool(13,ispec) + i14 = ibool(14,ispec) + i15 = ibool(15,ispec) + i16 = ibool(16,ispec) + i17 = ibool(17,ispec) + i18 = ibool(18,ispec) + i19 = ibool(19,ispec) + i20 = ibool(20,ispec) + i21 = ibool(21,ispec) + i22 = ibool(22,ispec) + i23 = ibool(23,ispec) + i24 = ibool(24,ispec) + i25 = ibool(25,ispec) + i26 = ibool(26,ispec) + i27 = ibool(27,ispec) + endif + + if (.true.) then + +! for the six faces below it is important to make sure we write the four points +! in an order for which the normal to the face points outwards + +! test face 1 (bottom) + if (z(i1) < limit .and. z(i2) < limit .and. z(i3) < limit .and. z(i4) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i4,i3,i2,i1 + else + write(24,"(10(i9,1x))") ispec,i4,i3,i2,i1,i11,i10,i9,i12,i21 + endif + endif + +! test face 2 (top) + if (z(i5) < limit .and. z(i6) < limit .and. z(i7) < limit .and. z(i8) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i5,i6,i7,i8 + else + write(24,"(10(i9,1x))") ispec,i5,i6,i7,i8,i17,i18,i19,i20,i26 + endif + endif + +! test face 3 (left) + if (z(i1) < limit .and. z(i4) < limit .and. z(i8) < limit .and. z(i5) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i5,i8,i4 + else + write(24,"(10(i9,1x))") ispec,i1,i5,i8,i4,i13,i20,i16,i12,i25 + endif + endif + +! test face 4 (right) + if (z(i2) < limit .and. z(i3) < limit .and. z(i7) < limit .and. z(i6) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i2,i3,i7,i6 + else + write(24,"(10(i9,1x))") ispec,i2,i3,i7,i6,i10,i15,i18,i14,i23 + endif + endif + +! test face 5 (front) + if (z(i1) < limit .and. z(i2) < limit .and. z(i6) < limit .and. z(i5) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i2,i6,i5 + else + write(24,"(10(i9,1x))") ispec,i1,i2,i6,i5,i9,i14,i17,i13,i22 + endif + endif + +! test face 6 (back) + if (z(i4) < limit .and. z(i3) < limit .and. z(i7) < limit .and. z(i8) < limit) then + if (NGNOD == 8) then + write(24,*) ispec,i3,i4,i8,i7 + else + write(24,"(10(i9,1x))") ispec,i3,i4,i8,i7,i11,i16,i19,i15,i24 + endif + endif + + endif + + enddo + + close(24) + + print *,'File "absorbing_surface_file_bottom" has been successfully created' + print * + +! ************* generate "free_or_absorbing_surface_file_zmax" ************* + +! first count the number of faces that are along that edge + + count_faces_found = 0 + +! Zmax + size_of_model = zmax - zmin + limit = zmax - SMALL_RELATIVE_VALUE*size_of_model + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + already_found_a_face = .false. + +! test face 1 (bottom) + if (z(i1) > limit .and. z(i2) > limit .and. z(i3) > limit .and. z(i4) > limit) then + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 2 (top) + if (z(i5) > limit .and. z(i6) > limit .and. z(i7) > limit .and. z(i8) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 3 (left) + if (z(i1) > limit .and. z(i4) > limit .and. z(i8) > limit .and. z(i5) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 4 (right) + if (z(i2) > limit .and. z(i3) > limit .and. z(i7) > limit .and. z(i6) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 5 (front) + if (z(i1) > limit .and. z(i2) > limit .and. z(i6) > limit .and. z(i5) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + +! test face 6 (back) + if (z(i4) > limit .and. z(i3) > limit .and. z(i7) > limit .and. z(i8) > limit) then + if (already_found_a_face) stop 'error: element with two faces on the same edge found!' + count_faces_found = count_faces_found + 1 + already_found_a_face = .true. + endif + + enddo + + print *,'found ',count_faces_found,' full faces on face Zmax' + +!----------------------------- + + open(unit=24,file='DATA_3D/free_or_absorbing_surface_file_zmax',status='unknown',action='write') + +! write the total number of face elements + write(24,*) count_faces_found + +! loop on the whole mesh + do ispec = 1,nspec + + i1 = ibool(1,ispec) + i2 = ibool(2,ispec) + i3 = ibool(3,ispec) + i4 = ibool(4,ispec) + i5 = ibool(5,ispec) + i6 = ibool(6,ispec) + i7 = ibool(7,ispec) + i8 = ibool(8,ispec) + + if (NGNOD == 27) then + i9 = ibool(9,ispec) + i10 = ibool(10,ispec) + i11 = ibool(11,ispec) + i12 = ibool(12,ispec) + i13 = ibool(13,ispec) + i14 = ibool(14,ispec) + i15 = ibool(15,ispec) + i16 = ibool(16,ispec) + i17 = ibool(17,ispec) + i18 = ibool(18,ispec) + i19 = ibool(19,ispec) + i20 = ibool(20,ispec) + i21 = ibool(21,ispec) + i22 = ibool(22,ispec) + i23 = ibool(23,ispec) + i24 = ibool(24,ispec) + i25 = ibool(25,ispec) + i26 = ibool(26,ispec) + i27 = ibool(27,ispec) + endif + +! for the six faces below it is important to make sure we write the four points +! in an order for which the normal to the face points outwards + +! test face 1 (bottom) + if (z(i1) > limit .and. z(i2) > limit .and. z(i3) > limit .and. z(i4) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i4,i3,i2,i1 + else + write(24,"(10(i9,1x))") ispec,i4,i3,i2,i1,i11,i10,i9,i12,i21 + endif + endif + +! test face 2 (top) + if (z(i5) > limit .and. z(i6) > limit .and. z(i7) > limit .and. z(i8) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i5,i6,i7,i8 + else + write(24,"(10(i9,1x))") ispec,i5,i6,i7,i8,i17,i18,i19,i20,i26 + endif + endif + +! test face 3 (left) + if (z(i1) > limit .and. z(i4) > limit .and. z(i8) > limit .and. z(i5) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i5,i8,i4 + else + write(24,"(10(i9,1x))") ispec,i1,i5,i8,i4,i13,i20,i16,i12,i25 + endif + endif + +! test face 4 (right) + if (z(i2) > limit .and. z(i3) > limit .and. z(i7) > limit .and. z(i6) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i2,i3,i7,i6 + else + write(24,"(10(i9,1x))") ispec,i2,i3,i7,i6,i10,i15,i18,i14,i23 + endif + endif + +! test face 5 (front) + if (z(i1) > limit .and. z(i2) > limit .and. z(i6) > limit .and. z(i5) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i1,i2,i6,i5 + else + write(24,"(10(i9,1x))") ispec,i1,i2,i6,i5,i9,i14,i17,i13,i22 + endif + endif + +! test face 6 (back) + if (z(i4) > limit .and. z(i3) > limit .and. z(i7) > limit .and. z(i8) > limit) then + if (NGNOD == 8) then + write(24,*) ispec,i3,i4,i8,i7 + else + write(24,"(10(i9,1x))") ispec,i3,i4,i8,i7,i11,i16,i19,i15,i24 + endif + endif + + enddo + + close(24) + + print *,'File "free_or_absorbing_surface_file_zmax" has been successfully created' + print * + + end program create_database_files_for_external_faces_of_the_3D_mesh + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_pressure_plane_wave_with_Hamming_window_00deg_for_crack.f90 b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_pressure_plane_wave_with_Hamming_window_00deg_for_crack.f90 new file mode 100644 index 000000000..6c62971d2 --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/create_pressure_plane_wave_with_Hamming_window_00deg_for_crack.f90 @@ -0,0 +1,47 @@ + + program create_Hamming_window + + implicit none + + integer, parameter :: NSOURCES = 1000 + + double precision, parameter :: xs_min = 0.020d0 + double precision, parameter :: xs_size = 0.010d0 + + double precision, parameter :: factor_max = 1.d10 + +! pi + double precision, parameter :: PI = 3.141592653589793d0 + + integer :: isource + + double precision :: x,hamming + + do isource = 1,NSOURCES + +! Hamming apodization window +! see e.g. http://docs.scipy.org/doc/numpy/reference/generated/numpy.hamming.html +! and http://www.mathworks.fr/fr/help/signal/ref/hamming.html + x = dble(isource - 1) / dble(NSOURCES - 1) + hamming = 0.54d0 - 0.46d0*cos(2*PI*x) + + write(*,*) '# source ',isource + write(*,*) 'source_surf = .false.' + write(*,*) 'xs = ',xs_min + x * xs_size + write(*,*) 'zs = 0.075' + write(*,*) 'source_type = 1' + write(*,*) 'time_function_type = 1' + write(*,*) 'name_of_source_file = YYYYYYYYYYYYYYYYYY' + write(*,*) 'burst_band_width = 0.' + write(*,*) 'f0 = 0.5d6' + write(*,*) 'tshift = 0.d0' + write(*,*) 'angleforce = 0.0' + write(*,*) 'Mxx = 1.' + write(*,*) 'Mzz = 1.' + write(*,*) 'Mxz = 0.' + write(*,*) 'factor = ',factor_max * hamming + + enddo + + end program create_Hamming_window + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_crack_water_20_mm_thick_et_maille_par_25_elements_donc_element_size_=_0.02_sur_25_=_0.00080_=_0.8_mm.txt b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_crack_water_20_mm_thick_et_maille_par_25_elements_donc_element_size_=_0.02_sur_25_=_0.00080_=_0.8_mm.txt new file mode 100644 index 000000000..0d5233a94 --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_crack_water_20_mm_thick_et_maille_par_25_elements_donc_element_size_=_0.02_sur_25_=_0.00080_=_0.8_mm.txt @@ -0,0 +1,11 @@ + +- cluster CPU + +- email SolTechnic + +========================= + +- maillage crack (water 20 mm thick et maille par 25 elements, donc element size = 0.02 sur 25 = 0.00080 = 0.8 mm) + +-relecture articles Sharham et Alexis + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_du_crack_newer_upwards_OK_avec_petit_mailleur_ecrit_par_Dimitri.pdf b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_du_crack_newer_upwards_OK_avec_petit_mailleur_ecrit_par_Dimitri.pdf new file mode 100644 index 000000000..9bfbd185c Binary files /dev/null and b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_du_crack_newer_upwards_OK_avec_petit_mailleur_ecrit_par_Dimitri.pdf differ diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_du_crack_older_downwards_OK_avec_petit_mailleur_ecrit_par_Dimitri.pdf b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_du_crack_older_downwards_OK_avec_petit_mailleur_ecrit_par_Dimitri.pdf new file mode 100644 index 000000000..c79601f37 Binary files /dev/null and b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/maillage_du_crack_older_downwards_OK_avec_petit_mailleur_ecrit_par_Dimitri.pdf differ diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/make_all_and_run_all_for_Y_shaped_crack_mesher.bash b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/make_all_and_run_all_for_Y_shaped_crack_mesher.bash new file mode 100644 index 000000000..970269d80 --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/make_all_and_run_all_for_Y_shaped_crack_mesher.bash @@ -0,0 +1,26 @@ +#!/bin/bash + +rm -f xtutu1 xtutu2 xtutu3 *__genmod.* + +ifort -std08 -implicitnone -warn all -check all -debug -g -O0 -fp-stack-check -traceback -ftrapuv -o xtutu1 mesh_a_Y_shaped_crack_oriented_upwards.f90 + +ifort -std08 -implicitnone -warn all -check all -debug -g -O0 -fp-stack-check -traceback -ftrapuv -o xtutu2 create_database_files_for_external_faces_of_the_3D_mesh.f90 + +ifort -std08 -implicitnone -warn all -check all -debug -g -O0 -fp-stack-check -traceback -ftrapuv -o xtutu3 create_OpenDX_files_to_view_the_3D_mesh.f90 + +rm -f *__genmod.* + +mkdir -p DATA +mkdir -p DATA_3D + +# run the two codes to create the 2D mesh, the 3D mesh, and the database files for SPECFEM3D_Cartesian +./xtutu1 +echo "creating the database files for the outer faces of the 3D..." +./xtutu2 +echo "Done." + +echo " " +echo "Now creating an OpenDX file to check the 3D mesh created..." +./xtutu3 +echo "Done." + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/mesh_a_Y_shaped_crack_oriented_downwards.f90 b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/mesh_a_Y_shaped_crack_oriented_downwards.f90 new file mode 100644 index 000000000..72c9533b3 --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/mesh_a_Y_shaped_crack_oriented_downwards.f90 @@ -0,0 +1,467 @@ + + program mesh_a_Y_shaped_crack + + implicit none + +! points defining the different regions to mesh + double precision, dimension(31) :: xnpgeo,znpgeo + + integer :: ip1,ip2,ip3,ip4,nx,nz,ispec,ier + + xnpgeo(01) = 0.d0 + znpgeo(01) = 0.d0 + + xnpgeo(02) = 18.6 + znpgeo(02) = 0.d0 + + xnpgeo(03) = 50 + znpgeo(03) = 0.d0 + + xnpgeo(04) = 0.d0 + znpgeo(04) = 15.d0 + + xnpgeo(05) = 18.6 + znpgeo(05) = 15.d0 + + xnpgeo(06) = 50 + znpgeo(06) = 15.d0 + + xnpgeo(07) = 0.d0 + znpgeo(07) = 19.87 + + xnpgeo(08) = 18.6 + znpgeo(08) = 19.87 + + xnpgeo(09) = 50 + znpgeo(09) = 19.87 + + xnpgeo(10) = 0 + znpgeo(10) = 25.37 + + xnpgeo(11) = 18.6 + znpgeo(11) = 25.37 + + xnpgeo(12) = 23.36 + znpgeo(12) = 23.53 + + xnpgeo(13) = 50 + znpgeo(13) = 23.53 + + xnpgeo(14) = 0 + znpgeo(14) = 29.46 + + xnpgeo(15) = 23.95 + znpgeo(15) = 29.46 + + xnpgeo(16) = 50 + znpgeo(16) = 29.46 + + xnpgeo(17) = 0 + znpgeo(17) = 35 + + xnpgeo(18) = 26.697 + znpgeo(18) = 35 + + xnpgeo(19) = 50 + znpgeo(19) = 35 + + xnpgeo(20) = 0 + znpgeo(20) = 43 + + xnpgeo(21) = 30.85 + znpgeo(21) = 43 + + xnpgeo(22) = 50 + znpgeo(22) = 43 + + xnpgeo(23) = 0 + znpgeo(23) = 50 + + xnpgeo(24) = 30.85 + znpgeo(24) = 50 + + xnpgeo(25) = 50 + znpgeo(25) = 50 + + xnpgeo(26) = 0 + znpgeo(26) = 65 + + xnpgeo(27) = 30.85 + znpgeo(27) = 65 + + xnpgeo(28) = 50 + znpgeo(28) = 65 + + xnpgeo(29) = 0 + znpgeo(29) = 85 + + xnpgeo(30) = 30.85 + znpgeo(30) = 85 + + xnpgeo(31) = 50 + znpgeo(31) = 85 + +! convert to millimiters + xnpgeo(:) = xnpgeo(:) / 1000.d0 + znpgeo(:) = znpgeo(:) / 1000.d0 + + ! user output + write(*,*) + write(*,*) 'Saving the mesh in Gnuplot format...' + write(*,*) + + open(unit=20,file='gridfile.gnu',status='unknown',iostat=ier) + if (ier /= 0 ) stop 'Error opening gnuplot file for writing: gridfile.gnu' + +! generate the mesh for all the regions + ispec = 0 + +! region 01 + ip1 = 1 + ip2 = 2 + ip3 = 5 + ip4 = 4 + nx = 32 + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 02 + ip1 = 2 + ip2 = 3 + ip3 = 6 + ip4 = 5 + nx = 32 + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 03 + ip1 = 4 + ip2 = 5 + ip3 = 8 + ip4 = 7 + nx = 32 + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 04 + ip1 = 5 + ip2 = 6 + ip3 = 9 + ip4 = 8 + nx = 32 + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 05 + ip1 = 7 + ip2 = 8 + ip3 = 11 + ip4 = 10 + nx = 32 + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 06 + ip1 = 8 + ip2 = 9 + ip3 = 13 + ip4 = 12 + nx = 32 + nz = 6 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 07 + ip1 = 8 + ip2 = 12 + ip3 = 15 + ip4 = 11 + nx = 6 + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 08 + ip1 = 10 + ip2 = 11 + ip3 = 15 + ip4 = 14 + nx = 32 + nz = 6 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 09 + ip1 = 12 + ip2 = 13 + ip3 = 16 + ip4 = 15 + nx = 32 + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 10 + ip1 = 14 + ip2 = 15 + ip3 = 18 + ip4 = 17 + nx = 32 + nz = 6 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 11 + ip1 = 15 + ip2 = 16 + ip3 = 19 + ip4 = 18 + nx = 32 + nz = 6 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 12 + ip1 = 17 + ip2 = 18 + ip3 = 21 + ip4 = 20 + nx = 32 + nz = 10 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 13 + ip1 = 18 + ip2 = 19 + ip3 = 22 + ip4 = 21 + nx = 32 + nz = 10 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 14 + ip1 = 20 + ip2 = 21 + ip3 = 24 + ip4 = 23 + nx = 32 + nz = 9 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 15 + ip1 = 21 + ip2 = 22 + ip3 = 25 + ip4 = 24 + nx = 32 + nz = 9 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 16 + ip1 = 23 + ip2 = 24 + ip3 = 27 + ip4 = 26 + nx = 32 + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 17 + ip1 = 24 + ip2 = 25 + ip3 = 28 + ip4 = 27 + nx = 32 + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 18 + ip1 = 26 + ip2 = 27 + ip3 = 30 + ip4 = 29 + nx = 32 + nz = 26 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + +! region 19 + ip1 = 27 + ip2 = 28 + ip3 = 31 + ip4 = 30 + nx = 32 + nz = 26 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + + print *,'total number of mesh elements generated = ',ispec + if (ispec /= (32 + 32) * (26*2 + 19*3) + 6*7) stop 'error in total number of mesh elements created' + + close(20) + + ! create a Gnuplot script to display the grid + open(unit=20,file='plot_gridfile.gnu',status='unknown',iostat=ier) + if (ier /= 0 ) stop 'Error saving plotgnu file' + + write(20,*) '#set term wxt' + write(20,*) + write(20,*) 'set term postscript portrait color solid "Helvetica" 22' + write(20,*) 'set output "maillage_du_crack_OK_avec_petit_mailleur_ecrit_par_Dimitri.ps"' + write(20,*) + + ! use same unit length on both X and Y axes + write(20,*) 'set size ratio -1' + + ! size of our model + write(20,*) 'set xrange [0:0.050]' + write(20,*) 'set yrange [0:0.085]' + + ! draw rectangles showing the water and steel layers + write(20,*) 'set object 1 rect from 0,0.065 to 0.050,0.085 fc rgb "#99FFFF" back' + write(20,*) 'set object 2 rect from 0,0.035 to 0.050,0.050 fc rgb "#888888" back' + write(20,*) 'set object 3 rect from 0,0 to 0.050,0.015 fc rgb "#888888" back' + + write(20,*) 'plot "gridfile.gnu" title "" w l lc black' + write(20,*) 'pause -1 "Hit any key..."' + close(20) + + write(*,*) 'Mesh saved in Gnuplot format...' + write(*,*) + + end program mesh_a_Y_shaped_crack + +!--------- + + subroutine generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,xnpgeo,znpgeo) + + implicit none + + integer :: ip1,ip2,ip3,ip4,nx,nz,ispec + +! points defining the different regions to mesh + double precision, dimension(31) :: xnpgeo,znpgeo + + double precision :: ratio_ix,ratio_iz,ratio_ixplus1,ratio_izplus1 + double precision :: x1,z1,x2,z2,x3,z3,x4,z4 + + integer :: ix,iz + + double precision, dimension(:,:), allocatable :: x,z + + allocate(x(0:nx,0:nz)) + allocate(z(0:nx,0:nz)) + + do iz = 0,nz-1 + do ix = 0,nx-1 + +! generate one more mesh element + ispec = ispec + 1 + + ratio_ix = ix / dble(nx) + ratio_iz = iz / dble(nz) + + ratio_ixplus1 = (ix+1) / dble(nx) + ratio_izplus1 = (iz+1) / dble(nz) + +! point 1 + call interpolate_bilinear(x1,z1,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,ratio_ix,ratio_iz) + +! point 2 + call interpolate_bilinear(x2,z2,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,ratio_ixplus1,ratio_iz) + +! point 3 + call interpolate_bilinear(x3,z3,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,ratio_ixplus1,ratio_izplus1) + +! point 4 + call interpolate_bilinear(x4,z4,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,ratio_ix,ratio_izplus1) + +! save the points created + x(ix,iz) = x1 + z(ix,iz) = z1 + + x(ix+1,iz) = x2 + z(ix+1,iz) = z2 + + x(ix+1,iz+1) = x3 + z(ix+1,iz+1) = z3 + + x(ix,iz+1) = x4 + z(ix,iz+1) = z4 + + enddo + enddo + + call save_gnuplot_file(nx,nz,x,z) + + deallocate(x) + deallocate(z) + + end subroutine generate_region + +!--------- + + subroutine interpolate_bilinear(x,z,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,gamma_interp_x,gamma_interp_z) + + implicit none + + integer :: ip1,ip2,ip3,ip4 + +! points defining the different regions to mesh + double precision, dimension(31) :: xnpgeo,znpgeo + + double precision :: gamma_interp_x,gamma_interp_z + double precision :: x,z + double precision :: val1,val2,val3,val4 + + ! interpolation rule + val1 = xnpgeo(ip1) + val2 = xnpgeo(ip2) + val3 = xnpgeo(ip3) + val4 = xnpgeo(ip4) + x = val1 * (1.d0-gamma_interp_x) * (1.d0-gamma_interp_z) + & + val2 * gamma_interp_x * (1.d0-gamma_interp_z) + & + val3 * gamma_interp_x * gamma_interp_z + & + val4 * (1.d0-gamma_interp_x) * gamma_interp_z + + val1 = znpgeo(ip1) + val2 = znpgeo(ip2) + val3 = znpgeo(ip3) + val4 = znpgeo(ip4) + z = val1 * (1.d0-gamma_interp_x) * (1.d0-gamma_interp_z) + & + val2 * gamma_interp_x * (1.d0-gamma_interp_z) + & + val3 * gamma_interp_x * gamma_interp_z + & + val4 * (1.d0-gamma_interp_x) * gamma_interp_z + + end subroutine interpolate_bilinear + +!------------------------- + +!!! this is src/meshfem2D/save_gnuplot_file.f90 from SPECFEM2D to display the mesh in Gnuplot format + + subroutine save_gnuplot_file(nx,nz,x,z) + +! creates a Gnuplot file that displays the grid + + implicit none + + integer :: nx,nz + double precision, dimension(0:nx,0:nz) :: x,z + + ! local parameters + integer :: ier,ili,icol + + ! draw horizontal lines of the grid + do ili=0,nz + do icol=0,nx-1 + write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili)) + write(20,*) sngl(x(icol+1,ili)),sngl(z(icol+1,ili)) + write(20,10) + enddo + enddo + + ! draw vertical lines of the grid + do icol=0,nx + do ili=0,nz-1 + write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili)) + write(20,*) sngl(x(icol,ili+1)),sngl(z(icol,ili+1)) + write(20,10) + enddo + enddo + +10 format('') + + end subroutine save_gnuplot_file + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/mesh_a_Y_shaped_crack_oriented_upwards.f90 b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/mesh_a_Y_shaped_crack_oriented_upwards.f90 new file mode 100644 index 000000000..7b4bcf6a1 --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/mesh_a_Y_shaped_crack_oriented_upwards.f90 @@ -0,0 +1,1120 @@ + + program mesh_a_Y_shaped_crack + + implicit none + +! spatial dimension of the mesh + integer, parameter :: NDIM = 2 + +! total number of spectral elements to create +!! DK DK half of the elements in the NX direction + integer, parameter :: nx_base = 48 ! 32 ! DK DK I have densified by 1.5 +!! DK DK number of elements in the vertical direction in the water layer + integer, parameter :: nz_water = 39 ! 26 ! DK DK I have densified by 1.5 + integer, parameter :: ngnod = 4 + integer, parameter :: nspec = (nx_base + nx_base) * (19 + 19 + 7 + 6 + 7 + 6 + 3 + 19 + nz_water) + 6*7 + integer, parameter :: NGLOB_MAX = nspec * ngnod + +! number of mesh elements in the third direction of the mesh +! (direction in which we extend the mesh identical to itself to convert it from 2D to 3D) + integer, parameter :: ny = nx_base + nx_base + double precision, parameter :: Ymin = 0.d0, Ymax = 0.050d0 + double precision, parameter :: deltaY = (Ymax - Ymin) / dble(ny) + +! points defining the different regions to mesh + double precision, dimension(31) :: xnpgeo,znpgeo + +! introduce split nodes on the crack or not (if not, the crack becomes invisible i.e. it does not exist any more +! and the layers that contain it become homogeneous and continuous; this is useful to provide a reference solution to compare to) + logical, parameter :: INTRODUCE_SPLIT_NODES_ON_THE_CRACK_CLEAN = .false. + + double precision, dimension(NGLOB_MAX) :: x,z + integer, dimension(nspec) :: imat + integer, dimension(ngnod,nspec) :: ibool + + integer :: ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,material_number,region_number,ipoin3D,ispec3D,iy +!!!!!!!!!!!! integer :: ispec_free_surface_found + integer :: ispec_absorbing_surface_found + double precision :: y + double precision :: geomtol + double precision :: big_X_offset + double precision :: angle_crack_to_horizontal,x1,z1,x2,z2,x3,z3 + +! additional arrays needed to sort the mesh points and remove the multiples + integer, dimension(:), allocatable :: locval,ind,ninseg,iglob,iwork + logical, dimension(:), allocatable :: ifseg + double precision, dimension(:), allocatable :: xp,yp,work + + integer :: nseg,ioff,iseg,ig,i,j,nglob,ia + integer :: ieoff,ilocnum + + double precision :: xmaxval,xminval,ymaxval,yminval,xtol,xtypdist + +! flags for Stacey absorbing boundaries + integer, parameter :: IBOTTOM = 1 + integer, parameter :: IRIGHT = 2 + integer, parameter :: ITOP = 3 + integer, parameter :: ILEFT = 4 + +! very large value + double precision, parameter :: HUGEVAL = 1.d+30 + +! relative mesh tolerance for fast global numbering + double precision, parameter :: SMALLVALTOL = 1.d-10 + + double precision, parameter :: PI = 3.141592653589793d0 + + xnpgeo(01) = 0.d0 + znpgeo(01) = 0.d0 + + xnpgeo(02) = 21.28 + znpgeo(02) = 0.d0 + + xnpgeo(03) = 50 + znpgeo(03) = 0.d0 + + xnpgeo(04) = 0.d0 + znpgeo(04) = 15.d0 + + xnpgeo(05) = 21.28 + znpgeo(05) = 15.d0 + + xnpgeo(06) = 50 + znpgeo(06) = 15.d0 + + xnpgeo(07) = 0.d0 + znpgeo(07) = 28.53 + + xnpgeo(08) = 21.28 + znpgeo(08) = 28.53 + + xnpgeo(09) = 50 + znpgeo(09) = 28.53 + + xnpgeo(10) = 0 + znpgeo(10) = 35. + + xnpgeo(11) = 24.75 + znpgeo(11) = 35. + + xnpgeo(12) = 50. + znpgeo(12) = 35. + + xnpgeo(13) = 0 + znpgeo(13) = 38.82 + + xnpgeo(14) = 27.09 + znpgeo(14) = 38.82 + + xnpgeo(15) = 50. + znpgeo(15) = 38.82 + + xnpgeo(16) = 0 + znpgeo(16) = 44.7 + + xnpgeo(17) = 27.7 + znpgeo(17) = 44.7 + + xnpgeo(18) = 31.87 + znpgeo(18) = 42.75 + + xnpgeo(19) = 50 + znpgeo(19) = 42.75 + + xnpgeo(20) = 0 + znpgeo(20) = 47.89 + + xnpgeo(21) = 30.35 + znpgeo(21) = 47.89 + + xnpgeo(22) = 50 + znpgeo(22) = 47.89 + + xnpgeo(23) = 0 + znpgeo(23) = 50 + + xnpgeo(24) = 30.35 + znpgeo(24) = 50 + + xnpgeo(25) = 50 + znpgeo(25) = 50 + + xnpgeo(26) = 0 + znpgeo(26) = 65 + + xnpgeo(27) = 30.35 + znpgeo(27) = 65 + + xnpgeo(28) = 50 + znpgeo(28) = 65 + + xnpgeo(29) = 0 + znpgeo(29) = 85 + + xnpgeo(30) = 30.35 + znpgeo(30) = 85 + + xnpgeo(31) = 50 + znpgeo(31) = 85 + +! convert to millimiters + xnpgeo(:) = xnpgeo(:) / 1000.d0 + znpgeo(:) = znpgeo(:) / 1000.d0 + +! generate the mesh for all the regions + ispec = 0 + ipoin = 0 + +! define a huge offset that is larger than the (Xmax - Xmin) horizontal size of the model +! in order to very easily create split nodes for the crack using a trick +! (the trick is that I shift the X coordinate of one side of the crack by a huge value, +! then perform the sorting of the mesh coordinates, and thus the two points are kept because one +! one is shifted from the other, and after sorting I shift the second point back to its original location). +! By doing so I automatically get split nodes for the two sides of the crack + if (INTRODUCE_SPLIT_NODES_ON_THE_CRACK_CLEAN) then + big_X_offset = 1000.d0 * (maxval(xnpgeo) - minval(xnpgeo)) + else + big_X_offset = 0.d0 + endif + + print * + x1 = xnpgeo(9) + z1 = znpgeo(9) + x2 = xnpgeo(14) + z2 = znpgeo(14) + x3 = xnpgeo(8) + z3 = znpgeo(8) + angle_crack_to_horizontal = acos(( (x1-x3)*(x2-x3) + (z1-z3)*(z2-z3) ) / & + ( sqrt((x1-x3)**2 + (z1-z3)**2) * sqrt((x2-x3)**2 + (z2-z3)**2) )) * 180.d0 / PI + print *,'The crack makes an angle of ',angle_crack_to_horizontal,' degrees with respect to the horizontal' + print *,'The crack makes an angle of ',90.d0 - angle_crack_to_horizontal,' degrees with respect to the vertical' + print * + + x1 = xnpgeo(8) + z1 = znpgeo(8) + x2 = (xnpgeo(17) + xnpgeo(18)) / 2.d0 + z2 = (znpgeo(17) + znpgeo(18)) / 2.d0 + print *,'The crack has an approximate total length of ',sqrt((x1-x2)**2 + (z1-z2)**2),' m' + print * + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!! this is the bottom layer (labeled Layer 4 in the InkScape figure) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + material_number = 5 + +! region 01 + region_number = 1 + ip1 = 1 + ip2 = 2 + ip3 = 5 + ip4 = 4 + nx = nx_base + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 02 + region_number = 2 + ip1 = 2 + ip2 = 3 + ip3 = 6 + ip4 = 5 + nx = nx_base + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!! this is Layer 3 in the InkScape figure +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + material_number = 4 + +! region 03 + region_number = 3 + ip1 = 4 + ip2 = 5 + ip3 = 8 + ip4 = 7 + nx = nx_base + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 04 + region_number = 4 + ip1 = 5 + ip2 = 6 + ip3 = 9 + ip4 = 8 + nx = nx_base + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 05 + region_number = 5 + ip1 = 7 + ip2 = 8 + ip3 = 11 + ip4 = 10 + nx = nx_base + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 06 + region_number = 6 + ip1 = 8 + ip2 = 9 + ip3 = 12 + ip4 = 11 + nx = nx_base + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!! this is Layer 2 in the InkScape figure +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + material_number = 3 + +! region 07 + region_number = 7 + ip1 = 10 + ip2 = 11 + ip3 = 14 + ip4 = 13 + nx = nx_base + nz = 6 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 08 + region_number = 8 + ip1 = 11 + ip2 = 12 + ip3 = 15 + ip4 = 14 + nx = nx_base + nz = 6 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 09 + region_number = 9 + ip1 = 13 + ip2 = 14 + ip3 = 17 + ip4 = 16 + nx = nx_base + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 10 + region_number = 10 + ip1 = 14 + ip2 = 18 + ip3 = 21 + ip4 = 17 + nx = 6 + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 11 + region_number = 11 + ip1 = 14 + ip2 = 15 + ip3 = 19 + ip4 = 18 + nx = nx_base + nz = 6 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 12 + region_number = 12 + ip1 = 16 + ip2 = 17 + ip3 = 21 + ip4 = 20 + nx = nx_base + nz = 6 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 13 + region_number = 13 + ip1 = 18 + ip2 = 19 + ip3 = 22 + ip4 = 21 + nx = nx_base + nz = 7 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 14 + region_number = 14 + ip1 = 20 + ip2 = 21 + ip3 = 24 + ip4 = 23 + nx = nx_base + nz = 3 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 15 + region_number = 15 + ip1 = 21 + ip2 = 22 + ip3 = 25 + ip4 = 24 + nx = nx_base + nz = 3 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!! this is Layer 1 in the InkScape figure +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + material_number = 2 + +! region 16 + region_number = 16 + ip1 = 23 + ip2 = 24 + ip3 = 27 + ip4 = 26 + nx = nx_base + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 17 + region_number = 17 + ip1 = 24 + ip2 = 25 + ip3 = 28 + ip4 = 27 + nx = nx_base + nz = 19 + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!! this is the top layer (the water layer) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + material_number = 1 + +! region 18 + region_number = 18 + ip1 = 26 + ip2 = 27 + ip3 = 30 + ip4 = 29 + nx = nx_base + nz = nz_water + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + +! region 19 + region_number = 19 + ip1 = 27 + ip2 = 28 + ip3 = 31 + ip4 = 30 + nx = nx_base + nz = nz_water + call generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + + print *,'total number of mesh elements generated = ',ispec + if (ispec /= nspec) stop 'error in total number of mesh elements created' + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!---- create global mesh numbering, removing the multiples + + allocate(locval(NGLOB_MAX)) + allocate(ind(NGLOB_MAX)) + allocate(ninseg(NGLOB_MAX)) + allocate(iglob(NGLOB_MAX)) + allocate(ifseg(NGLOB_MAX)) + allocate(xp(NGLOB_MAX)) + allocate(yp(NGLOB_MAX)) + allocate(work(NGLOB_MAX)) + allocate(iwork(NGLOB_MAX)) + +! copy coordinates of the grid points (since sorting will destroy the array order, and we will need to access the original) + xp(:) = x(:) + yp(:) = z(:) + +! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +! Establish initial pointers + do ispec = 1,nspec + ieoff = ngnod*(ispec -1) + do ia = 1,ngnod + locval (ia+ieoff) = ia+ieoff + enddo + enddo + +! set up a local geometric tolerance + + xtypdist = +HUGEVAL + + do ispec = 1,nspec + + xminval = +HUGEVAL + yminval = +HUGEVAL + xmaxval = -HUGEVAL + ymaxval = -HUGEVAL + ieoff = ngnod*(ispec-1) + do ilocnum = 1,ngnod + xmaxval = max(xp(ieoff+ilocnum),xmaxval) + xminval = min(xp(ieoff+ilocnum),xminval) + ymaxval = max(yp(ieoff+ilocnum),ymaxval) + yminval = min(yp(ieoff+ilocnum),yminval) + enddo + +! compute the minimum typical "size" of an element in the mesh + xtypdist = min(xtypdist,xmaxval-xminval) + xtypdist = min(xtypdist,ymaxval-yminval) + + enddo + +! define a tolerance, small with respect to the minimum size + xtol = SMALLVALTOL * xtypdist + +! print *,'DK DK xtol = ',xtol + + ifseg(:) = .false. + nseg = 1 + ifseg(1) = .true. + ninseg(1) = NGLOB_MAX + + do j = 1,NDIM +! Sort within each segment + ioff = 1 + do iseg = 1,nseg + if (j == 1) then + call rank (xp(ioff),ind,ninseg(iseg)) + else + call rank (yp(ioff),ind,ninseg(iseg)) + endif + call swap(xp(ioff),work,ind,ninseg(iseg)) + call swap(yp(ioff),work,ind,ninseg(iseg)) + call iswap(locval(ioff),iwork,ind,ninseg(iseg)) + ioff = ioff + ninseg(iseg) + enddo +! Check for jumps in current coordinate + if (j == 1) then + do i = 2,NGLOB_MAX + if (abs(xp(i)-xp(i-1)) > xtol) ifseg(i) = .true. + enddo + else + do i = 2,NGLOB_MAX + if (abs(yp(i)-yp(i-1)) > xtol) ifseg(i) = .true. + enddo + endif +! Count up number of different segments + nseg = 0 + do i = 1,NGLOB_MAX + if (ifseg(i)) then + nseg = nseg + 1 + ninseg(nseg) = 1 + else + ninseg(nseg) = ninseg(nseg) + 1 + endif + enddo + enddo +! +! Assign global node numbers (now sorted lexicographically!) +! + ig = 0 + do i = 1,NGLOB_MAX + if (ifseg(i)) then + ig = ig+1 +!! DK DK sep 2018: added this to save the list of grid points without multiples + x(ig) = xp(i) + z(ig) = yp(i) +!! DK DK sep 2018: added this to save the list of grid points without multiples + +!! DK DK restore the correct physical location of the split nodes, i.e. remove the shift in X that I have introduced +!! DK DK purposely in order to create split nodes automatically in the mesh point sorting routine and point multiple removal + if (INTRODUCE_SPLIT_NODES_ON_THE_CRACK_CLEAN) then + if (x(ig) < -SMALLVALTOL) x(ig) = x(ig) + big_X_offset + if (x(ig) > big_X_offset) x(ig) = x(ig) - big_X_offset + endif + + endif + iglob(locval(i)) = ig + enddo + + nglob = ig + +! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +! get result in my format + do ispec = 1,nspec + ieoff = ngnod*(ispec - 1) + ilocnum = 0 + do ia = 1,ngnod + ilocnum = ilocnum + 1 + ibool(ia,ispec) = iglob(ilocnum + ieoff) + enddo + enddo + + deallocate(locval) + deallocate(ind) + deallocate(ninseg) + deallocate(iglob) + deallocate(ifseg) + deallocate(xp) + deallocate(yp) + deallocate(work) + deallocate(iwork) + +! check the numbering obtained + if (minval(ibool) /= 1 .or. maxval(ibool) /= nglob) stop 'Error while generating global numbering' + + print *,'total number of unique mesh points generated = ',nglob + if (nglob > NGLOB_MAX) stop 'error in total number of unique mesh points created' + +! save the mesh files in the format needed by SPECFEM2D for external meshes + open(unit=20,file='DATA/mesh_file',status='unknown') + write(20,*) nspec + do ispec = 1,nspec + write(20,*) (ibool(ia,ispec),ia=1,ngnod) + enddo + close(20) + + open(unit=20,file='DATA/nodes_coords_file',status='unknown') + write(20,*) nglob + do ipoin = 1,nglob + write(20,*) x(ipoin),z(ipoin) + enddo + close(20) + + open(unit=20,file='DATA/materials_file',status='unknown') + do ispec = 1,nspec + write(20,*) imat(ispec) + enddo + close(20) + +! geometrical tolerance to detect the edges of the model (1/1000th of the total model size) + geomtol = 0.085d0 / 1000.d0 + + open(unit=20,file='DATA/free_surface_file',status='unknown') + write(20,*) '0' ! create an empty file (file listing zero elements) +! write(20,*) nx_base + nx_base ! elements at the top free surface +! ispec_free_surface_found = 0 +! do ispec = 1,nspec +! check if the element is located in the water +! and if the top edge of this element belongs to the acoustic free surface, which is located at Z = 0.085 +! if (imat(ispec) == 1 .and. z(ibool(3,ispec)) > 0.085d0 - geomtol .and. z(ibool(4,ispec)) > 0.085d0 - geomtol) then +! ispec_free_surface_found = ispec_free_surface_found + 1 +! ! 'acoustic_surface' contains 1/ element number, 2/ number of nodes that form the free surface, +! ! 3/ first node on the free surface, 4/ second node on the free surface, if relevant (if 2/ is equal to 2) +! write(20,*) ispec,' 2 ',ibool(3,ispec),ibool(4,ispec) +! endif +! enddo +! if (ispec_free_surface_found /= nx_base + nx_base) stop 'error in detecting acoustic free surface' + close(20) + + open(unit=20,file='DATA/axial_elements_file',status='unknown') + write(20,*) '0' ! create an empty file (file listing zero elements) + close(20) + + open(unit=20,file='DATA/absorbing_surface_file',status='unknown') +! write(20,*) '0' ! create an empty file (file listing zero elements) + +! count the number of absorbing edges + ispec_absorbing_surface_found = 0 + do ispec = 1,nspec + if (x(ibool(1,ispec)) < geomtol .and. x(ibool(4,ispec)) < geomtol) & + ispec_absorbing_surface_found = ispec_absorbing_surface_found + 1 + if (x(ibool(2,ispec)) > 0.050d0 - geomtol .and. x(ibool(3,ispec)) > 0.050d0 - geomtol) & + ispec_absorbing_surface_found = ispec_absorbing_surface_found + 1 + + if (z(ibool(1,ispec)) < geomtol .and. z(ibool(2,ispec)) < geomtol) & + ispec_absorbing_surface_found = ispec_absorbing_surface_found + 1 + if (z(ibool(3,ispec)) > 0.085d0 - geomtol .and. z(ibool(4,ispec)) > 0.085d0 - geomtol) & + ispec_absorbing_surface_found = ispec_absorbing_surface_found + 1 + enddo + write(20,*) ispec_absorbing_surface_found + +! write the absorbing edges +! ! 'absorbing_surface' contains 1/ element number, 2/ number of nodes that form the free surface, +! ! 3/ first node on the free surface, 4/ second node on the free surface, if relevant (if 2/ is equal to 2) + do ispec = 1,nspec + if (x(ibool(1,ispec)) < geomtol .and. x(ibool(4,ispec)) < geomtol) then + write(20,*) ispec,' 2 ',ibool(1,ispec),ibool(4,ispec),ILEFT + endif + + if (x(ibool(2,ispec)) > 0.050d0 - geomtol .and. x(ibool(3,ispec)) > 0.050d0 - geomtol) then + write(20,*) ispec,' 2 ',ibool(2,ispec),ibool(3,ispec),IRIGHT + endif + + if (z(ibool(1,ispec)) < geomtol .and. z(ibool(2,ispec)) < geomtol) then + write(20,*) ispec,' 2 ',ibool(1,ispec),ibool(2,ispec),IBOTTOM + endif + + if (z(ibool(3,ispec)) > 0.085d0 - geomtol .and. z(ibool(4,ispec)) > 0.085d0 - geomtol) then + write(20,*) ispec,' 2 ',ibool(3,ispec),ibool(4,ispec),ITOP + endif + enddo + + close(20) + + open(unit=20,file='DATA/Surf_acforcing_Bottom_enforcing_mesh',status='unknown') + write(20,*) '0' ! create an empty file (file listing zero elements) + close(20) + + open(unit=20,file='DATA/absorbing_cpml_file',status='unknown') + write(20,*) '0' ! create an empty file (file listing zero elements) + close(20) + + open(unit=20,file='DATA/courbe_eros_nodes',status='unknown') + write(20,*) '0' ! create an empty file (file listing zero elements) + close(20) + + +! save the mesh in Gnuplot format for verification + call save_gnuplot_file(x,z,nglob,ibool,ngnod,nspec) + + print *,'done creating and saving the 2D mesh' + +!! DK DK +!! DK DK also save a 3D mesh, by extending the mesh identical to itself in the third direction +!! DK DK + + print *,'creating and saving the 3D mesh...' + print *,'the 3D mesh contains ',nglob * (ny + 1),' unique points, and ',nspec * ny,' elements' + +! create SPECFEM3D_Cartesian mesh file for the points + open(unit=20,file='DATA_3D/nodes_coords_file',status='unknown',action='write') + write(20,*) nglob * (ny + 1) + ipoin3D = 0 + do iy = 1,ny + 1 + y = (iy - 1) * deltaY + do ipoin = 1,nglob + ipoin3D = ipoin3D + 1 + write(20,*) ipoin3D,x(ipoin),y,z(ipoin) + enddo + enddo + close(20) + + open(unit=20,file='DATA_3D/mesh_file',status='unknown') + write(20,*) nspec * ny + ispec3D = 0 + do iy = 1,ny + do ispec = 1,nspec + ispec3D = ispec3D + 1 +! the 3D element is created by using the 2D element, and the 2D element of the next plane in the Y direction, +! which is nglob grid points further in the list of points, due to the fact that to create the 3D mesh +! we extend of the 2D mesh (which contains nglob points) identical to itself in the third direction. +! To get the point order right, see file numbering_convention_8_nodes.png + write(20,*) ispec3D,ibool(1,ispec) + (iy-1) * nglob,ibool(2,ispec) + (iy-1) * nglob, & + ibool(2,ispec) + iy * nglob,ibool(1,ispec) + iy * nglob, & + ibool(4,ispec) + (iy-1) * nglob,ibool(3,ispec) + (iy-1) * nglob, & + ibool(3,ispec) + iy * nglob,ibool(4,ispec) + iy * nglob + enddo + enddo + close(20) + + open(unit=20,file='DATA_3D/materials_file',status='unknown') + ispec3D = 0 + do iy = 1,ny + do ispec = 1,nspec + ispec3D = ispec3D + 1 + write(20,*) ispec3D,imat(ispec) + enddo + enddo + close(20) + + print *,'done creating and saving the 3D mesh' + + end program mesh_a_Y_shaped_crack + +!--------- + + subroutine generate_region(ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,xnpgeo,znpgeo,material_number,x,z,imat,nspec,NGLOB_MAX, & + region_number,big_X_offset) + + implicit none + + integer :: ip1,ip2,ip3,ip4,nx,nz,ispec,ipoin,material_number,region_number + double precision :: big_X_offset + +! points defining the different regions to mesh + double precision, dimension(31) :: xnpgeo,znpgeo + + double precision :: ratio_ix,ratio_iz,ratio_ixplus1,ratio_izplus1 + double precision :: x1,z1,x2,z2,x3,z3,x4,z4 + double precision :: xbackup1,zbackup1 + double precision :: xbackup2 + double precision :: big_offset_local + + integer :: ix,iz + + integer :: nspec,NGLOB_MAX + double precision, dimension(NGLOB_MAX) :: x,z + integer, dimension(nspec) :: imat + +! introduce split nodes on the crack or not (if not, the crack becomes invisible i.e. it does not exist any more +! and the layers that contain it become homogeneous and continuous; this is useful to provide a reference solution to compare to) + logical, parameter :: INTRODUCE_SPLIT_NODES_ON_THE_CRACK_UGLY = .true. + + if (INTRODUCE_SPLIT_NODES_ON_THE_CRACK_UGLY) then +! big_offset_local = 0.03d0 * (50.d0/1000.d0/64.d0) !! 3% of average size of an element + big_offset_local = 0.0003d0 * (50.d0/1000.d0/64.d0) !! 0.03% of average size of an element + else + big_offset_local = 0.d0 + endif + +!! DK DK YYYYYYYYYYYYYYYYYYYYYY + if (INTRODUCE_SPLIT_NODES_ON_THE_CRACK_UGLY) then + if (region_number == 5) then + xbackup1 = xnpgeo(ip3) + xnpgeo(ip3) = xnpgeo(ip3) - big_offset_local + endif + + if (region_number == 7) then + xbackup1 = xnpgeo(ip2) + xnpgeo(ip2) = xnpgeo(ip2) - big_offset_local + xbackup2 = xnpgeo(ip3) + xnpgeo(ip3) = xnpgeo(ip3) - big_offset_local + endif + + if (region_number == 9) then + xbackup1 = xnpgeo(ip2) + xnpgeo(ip2) = xnpgeo(ip2) - big_offset_local + endif + + if (region_number == 10) then + xbackup1 = xnpgeo(ip1) + zbackup1 = znpgeo(ip1) + xnpgeo(ip1) = xnpgeo(ip1) + big_offset_local + znpgeo(ip1) = znpgeo(ip1) + big_offset_local + endif + endif +!! DK DK YYYYYYYYYYYYYYYYYYYYYY + + do iz = 0,nz-1 + do ix = 0,nx-1 + +! generate one more mesh element + ispec = ispec + 1 + +! store the material number + imat(ispec) = material_number + + ratio_ix = ix / dble(nx) + ratio_iz = iz / dble(nz) + + ratio_ixplus1 = (ix+1) / dble(nx) + ratio_izplus1 = (iz+1) / dble(nz) + +! point 1 + call interpolate_bilinear(x1,z1,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,ratio_ix,ratio_iz) + +! point 2 + call interpolate_bilinear(x2,z2,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,ratio_ixplus1,ratio_iz) + +! point 3 + call interpolate_bilinear(x3,z3,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,ratio_ixplus1,ratio_izplus1) + +! point 4 + call interpolate_bilinear(x4,z4,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,ratio_ix,ratio_izplus1) + +! shift the points to introduce split nodes automatically if needed + if (region_number == 5 .or. region_number == 7 .or. region_number == 9) then + ! right edge of the region left of the crack, shift the two points that are on the IRIGHT edge of the element, + ! which corresponds to the crack + if (ix == nx-1) then + x2 = x2 + big_X_offset + x3 = x3 + big_X_offset + endif + endif + + if (region_number == 11) then + ! left edge of the region right of the crack, shift the two points that are on the ILEFT edge of the element, + ! which corresponds to the crack + if (ix == 0) then + ! use a negative sign this time, because the center point of the Y of the crack needs to be split twice, + ! once to the right and the second time to the left, since that particular point is shared between three branches + ! of the crack instead of shared between only two side + x1 = x1 - big_X_offset + x4 = x4 - big_X_offset + endif + endif + +! save the points created + x(ipoin+1) = x1 + z(ipoin+1) = z1 + + x(ipoin+2) = x2 + z(ipoin+2) = z2 + + x(ipoin+3) = x3 + z(ipoin+3) = z3 + + x(ipoin+4) = x4 + z(ipoin+4) = z4 + + ipoin = ipoin + 4 + + enddo + enddo + +!! DK DK YYYYYYYYYYYYYYYYYYYYYY +! restore the right geometrical position of the nodes, after having introduced the splitting above + if (INTRODUCE_SPLIT_NODES_ON_THE_CRACK_UGLY) then + if (region_number == 5) then + xnpgeo(ip3) = xbackup1 + endif + + if (region_number == 7) then + xnpgeo(ip2) = xbackup1 + xnpgeo(ip3) = xbackup2 + endif + + if (region_number == 9) then + xnpgeo(ip2) = xbackup1 + endif + + if (region_number == 10) then + xnpgeo(ip1) = xbackup1 + znpgeo(ip1) = zbackup1 + endif + endif +!! DK DK YYYYYYYYYYYYYYYYYYYYYY + + end subroutine generate_region + +!--------- + + subroutine interpolate_bilinear(x,z,ip1,ip2,ip3,ip4,xnpgeo,znpgeo,gamma_interp_x,gamma_interp_z) + + implicit none + + integer :: ip1,ip2,ip3,ip4 + +! points defining the different regions to mesh + double precision, dimension(31) :: xnpgeo,znpgeo + + double precision :: gamma_interp_x,gamma_interp_z + double precision :: x,z + double precision :: val1,val2,val3,val4 + + ! interpolation rule + val1 = xnpgeo(ip1) + val2 = xnpgeo(ip2) + val3 = xnpgeo(ip3) + val4 = xnpgeo(ip4) + x = val1 * (1.d0-gamma_interp_x) * (1.d0-gamma_interp_z) + & + val2 * gamma_interp_x * (1.d0-gamma_interp_z) + & + val3 * gamma_interp_x * gamma_interp_z + & + val4 * (1.d0-gamma_interp_x) * gamma_interp_z + + val1 = znpgeo(ip1) + val2 = znpgeo(ip2) + val3 = znpgeo(ip3) + val4 = znpgeo(ip4) + z = val1 * (1.d0-gamma_interp_x) * (1.d0-gamma_interp_z) + & + val2 * gamma_interp_x * (1.d0-gamma_interp_z) + & + val3 * gamma_interp_x * gamma_interp_z + & + val4 * (1.d0-gamma_interp_x) * gamma_interp_z + + end subroutine interpolate_bilinear + +!------------------------- + +!!! this is adapted from src/meshfem2D/save_gnuplot_file.f90 of SPECFEM2D to display the mesh in Gnuplot format + + subroutine save_gnuplot_file(x,z,nglob,ibool,ngnod,nspec) + +! creates a Gnuplot file that displays the grid + + implicit none + + integer :: ngnod,nspec,nglob + integer, dimension(ngnod,nspec) :: ibool + double precision, dimension(nglob) :: x,z + + ! local parameters + integer :: ispec,ier + + open(unit=20,file='gridfile.gnu',status='unknown',iostat=ier) + if (ier /= 0) stop 'Error opening gnuplot file for writing: gridfile.gnu' + + do ispec = 1,nspec + + ! draw horizontal lines of the grid + write(20,*) sngl(x(ibool(1,ispec))),sngl(z(ibool(1,ispec))) + write(20,*) sngl(x(ibool(2,ispec))),sngl(z(ibool(2,ispec))) + write(20,10) + + write(20,*) sngl(x(ibool(3,ispec))),sngl(z(ibool(3,ispec))) + write(20,*) sngl(x(ibool(4,ispec))),sngl(z(ibool(4,ispec))) + write(20,10) + + ! draw vertical lines of the grid + write(20,*) sngl(x(ibool(1,ispec))),sngl(z(ibool(1,ispec))) + write(20,*) sngl(x(ibool(4,ispec))),sngl(z(ibool(4,ispec))) + write(20,10) + + write(20,*) sngl(x(ibool(2,ispec))),sngl(z(ibool(2,ispec))) + write(20,*) sngl(x(ibool(3,ispec))),sngl(z(ibool(3,ispec))) + write(20,10) + + enddo + + close(20) + + 10 format('') + + ! create a Gnuplot script to display the grid + open(unit=20,file='plot_gridfile.gnu',status='unknown',iostat=ier) + if (ier /= 0) stop 'Error saving plot_gridfile.gnu file' + + write(20,*) '#set term wxt' + write(20,*) 'set term x11' + write(20,*) 'set term qt' + write(20,*) + write(20,*) '#set term postscript portrait color solid "Helvetica" 22' + write(20,100) '#set output "maillage_du_crack_newer_upwards_OK_avec_petit_mailleur_ecrit_par_Dimitri.ps"' + write(20,*) + + ! use same unit length on both X and Y axes + write(20,*) 'set size ratio -1' + + ! size of our model + write(20,*) 'set xrange [0:0.050]' + write(20,*) 'set yrange [0:0.085]' + + ! draw rectangles showing the water and steel layers + write(20,*) 'set object 1 rect from 0,0.065 to 0.050,0.085 fc rgb "#99FFFF" back' + write(20,*) 'set object 2 rect from 0,0.035 to 0.050,0.050 fc rgb "#888888" back' + write(20,*) 'set object 3 rect from 0,0 to 0.050,0.015 fc rgb "#888888" back' + + write(20,*) 'plot "gridfile.gnu" title "" w l lc black' + write(20,*) 'pause -1 "Hit any key..."' + close(20) + + 100 format(a200) + + end subroutine save_gnuplot_file + +! +!----------------------------------------------------------------------- +! + + subroutine rank(A,IND,N) +! +! Use Heap Sort (p 233 Numerical Recipes) +! + implicit none + + integer N + double precision A(N) + integer IND(N) + + integer i,j,l,ir,indx + double precision q + + do j = 1,N + IND(j)=j + enddo + + if (n == 1) return + L=n/2+1 + ir=n + 100 continue + if (l > 1) then + l=l-1 + indx=ind(l) + q=a(indx) + ELSE + indx=ind(ir) + q=a(indx) + ind(ir)=ind(1) + ir=ir-1 + if (ir == 1) then + ind(1)=indx + return + endif + endif + i=l + j=l+l + 200 continue + if (J <= IR) then + if (J < IR) then + if (A(IND(j)) < A(IND(j+1))) j=j+1 + endif + if (q < A(IND(j))) then + IND(I)=IND(J) + I=J + J=J+J + ELSE + J=IR+1 + endif + goto 200 + endif + IND(I)=INDX + goto 100 + + end subroutine rank + +!----------------------------------------------------------------------- + + subroutine swap(a,w,ind,n) +! +! Use IND to sort array A (p 233 Numerical Recipes) +! + implicit none + + integer n + double precision A(N),W(N) + integer IND(N) + + integer j + + W(:) = A(:) + + do j = 1,N + A(j) = W(ind(j)) + enddo + + end subroutine swap + +!----------------------------------------------------------------------- + + subroutine iswap(a,w,ind,n) +! +! Use IND to sort array A +! + implicit none + + integer n + integer A(N),W(N),IND(N) + + integer j + + W(:) = A(:) + + do j = 1,N + A(j) = W(ind(j)) + enddo + + end subroutine iswap + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_downwards.svg b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_downwards.svg new file mode 100644 index 000000000..fca63072e --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_downwards.svg @@ -0,0 +1,496 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + y + z + 1 + + 2 + + 3 + + 4 + + Water + + + + + + + + + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_downwards_with_19_regions.png b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_downwards_with_19_regions.png new file mode 100644 index 000000000..8842b9bc9 Binary files /dev/null and b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_downwards_with_19_regions.png differ diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards.svg b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards.svg new file mode 100644 index 000000000..e40df799c --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards.svg @@ -0,0 +1,496 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + y + z + 1 + + 2 + + 3 + + 4 + + Water + + + + + + + + + diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_19_regions_better.png b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_19_regions_better.png new file mode 100644 index 000000000..821a7aa28 Binary files /dev/null and b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_19_regions_better.png differ diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_19_regions_even_better_final_chosen.png b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_19_regions_even_better_final_chosen.png new file mode 100644 index 000000000..b151804c3 Binary files /dev/null and b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_19_regions_even_better_final_chosen.png differ diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_21_regions_a_bit_too_complex.png b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_21_regions_a_bit_too_complex.png new file mode 100644 index 000000000..03d7f186c Binary files /dev/null and b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/setting_2D_crack_to_extend_to_3D_oriented_upwards_with_21_regions_a_bit_too_complex.png differ diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/visualize_mesh_DX.cfg b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/visualize_mesh_DX.cfg new file mode 100644 index 000000000..26afbdf71 --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/visualize_mesh_DX.cfg @@ -0,0 +1,45 @@ +// +// time: Fri Sep 21 02:12:54 2018 +// +// version: 3.2.0 (format), 4.4.4 (DX) +// +// +// panel[0]: position = (0.0750,0.6981), size = 0.2781x0.1611, startup = 1, devstyle = 1 +// title: value = Control Panel +// +// workspace: width = 500, height = 500 +// layout: snap = 0, width = 50, height = 50, align = NN +// +// panel[2]: position = (0.3609,0.6981), size = 0.2776x0.1611, startup = 1, devstyle = 1 +// title: value = Control Panel +// +// workspace: width = 500, height = 500 +// layout: snap = 0, width = 50, height = 50, align = NN +// +// node Image[2]: +// depth: value = 24 +// window: position = (0.1562,0.0426), size = 0.6641x0.8815 +// input[1]: defaulting = 0, value = "Image_2" +// input[4]: defaulting = 0, value = 1 +// input[5]: defaulting = 0, value = [0.025 0.025 0.0425] +// input[6]: defaulting = 0, value = [0.0943607 -0.252846 0.0756326] +// input[7]: defaulting = 0, value = 0.154491 +// input[8]: defaulting = 0, value = 1261 +// input[9]: defaulting = 0, value = 0.722839 +// input[10]: defaulting = 0, value = [0.0387653 0.127856 0.991035] +// input[11]: defaulting = 1, value = 30.0001 +// input[12]: defaulting = 0, value = 0 +// input[14]: defaulting = 0, value = 1 +// input[15]: defaulting = 0, value = "none" +// input[16]: defaulting = 0, value = "none" +// input[17]: defaulting = 0, value = 1 +// input[18]: defaulting = 0, value = 1 +// input[19]: defaulting = 0, value = 1 +// input[22]: defaulting = 0, value = "black" +// input[25]: defaulting = 0, value = "/home/latychev/boo" +// input[26]: defaulting = 0, value = "eps gray dpi=86 orient=portrait" +// input[29]: defaulting = 0, value = 0 +// input[33]: defaulting = 0, value = 1 +// input[36]: defaulting = 0, value = 1 +// input[41]: defaulting = 0, value = "rotate" +// internal caching: 1 diff --git a/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/visualize_mesh_DX.net b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/visualize_mesh_DX.net new file mode 100644 index 000000000..9cd10ae65 --- /dev/null +++ b/utils/mesher_analytical_Dimitri_for_2D_and_3D_cracks/visualize_mesh_DX.net @@ -0,0 +1,538 @@ +// +// time: Fri Sep 21 02:12:54 2018 +// +// version: 3.2.0 (format), 4.4.4 (DX) +// +// +// MODULE main +// workspace: width = 393, height = 599 +// layout: snap = 0, width = 50, height = 50, align = NN +// +macro main( +) -> ( +) { + // + // node AmbientLight[1]: x = 303, y = 277, inputs = 1, label = AmbientLight + // +main_AmbientLight_1_out_1 = + AmbientLight( + main_AmbientLight_1_in_1 + ) [instance: 1, cache: 1]; + // + // node Import[2]: x = 193, y = 42, inputs = 6, label = Import + // input[1]: defaulting = 0, visible = 1, type = 32, value = "DX_fullmesh.dx" + // input[3]: defaulting = 0, visible = 1, type = 32, value = "dx" + // +main_Import_2_out_1 = + Import( + main_Import_2_in_1, + main_Import_2_in_2, + main_Import_2_in_3, + main_Import_2_in_4, + main_Import_2_in_5, + main_Import_2_in_6 + ) [instance: 2, cache: 1]; + // + // node AutoColor[1]: x = 178, y = 167, inputs = 10, label = AutoColor + // +main_AutoColor_1_out_1, +main_AutoColor_1_out_2 = + AutoColor( + main_Import_2_out_1, + main_AutoColor_1_in_2, + main_AutoColor_1_in_3, + main_AutoColor_1_in_4, + main_AutoColor_1_in_5, + main_AutoColor_1_in_6, + main_AutoColor_1_in_7, + main_AutoColor_1_in_8, + main_AutoColor_1_in_9, + main_AutoColor_1_in_10 + ) [instance: 1, cache: 1]; + // + // node ShowConnections[1]: x = 49, y = 315, inputs = 1, label = ShowConnections + // +main_ShowConnections_1_out_1 = + ShowConnections( + main_AutoColor_1_out_1 + ) [instance: 1, cache: 1]; + // + // node Collect[1]: x = 167, y = 429, inputs = 4, label = Collect + // +main_Collect_1_out_1 = + Collect( + main_ShowConnections_1_out_1, + main_AutoColor_1_out_1, + main_AmbientLight_1_out_1, + main_Collect_1_in_4 + ) [instance: 1, cache: 1]; + // + // node Image[2]: x = 192, y = 540, inputs = 49, label = Image + // input[1]: defaulting = 0, visible = 0, type = 67108863, value = "Image_2" + // input[4]: defaulting = 0, visible = 0, type = 1, value = 1 + // input[5]: defaulting = 0, visible = 0, type = 8, value = [0.025 0.025 0.0425] + // input[6]: defaulting = 0, visible = 0, type = 8, value = [0.0943607 -0.252846 0.0756326] + // input[7]: defaulting = 0, visible = 0, type = 5, value = 0.154491 + // input[8]: defaulting = 0, visible = 0, type = 1, value = 1261 + // input[9]: defaulting = 0, visible = 0, type = 5, value = 0.722839 + // input[10]: defaulting = 0, visible = 0, type = 8, value = [0.0387653 0.127856 0.991035] + // input[11]: defaulting = 1, visible = 0, type = 5, value = 30.0001 + // input[12]: defaulting = 0, visible = 0, type = 1, value = 0 + // input[14]: defaulting = 0, visible = 0, type = 1, value = 1 + // input[15]: defaulting = 0, visible = 0, type = 32, value = "none" + // input[16]: defaulting = 0, visible = 0, type = 32, value = "none" + // input[17]: defaulting = 0, visible = 0, type = 1, value = 1 + // input[18]: defaulting = 0, visible = 0, type = 1, value = 1 + // input[19]: defaulting = 0, visible = 0, type = 1, value = 1 + // input[22]: defaulting = 0, visible = 0, type = 32, value = "black" + // input[25]: defaulting = 0, visible = 0, type = 32, value = "/home/latychev/boo" + // input[26]: defaulting = 0, visible = 0, type = 32, value = "eps gray dpi=86 orient=portrait" + // input[29]: defaulting = 0, visible = 0, type = 3, value = 0 + // input[33]: defaulting = 0, visible = 0, type = 3, value = 1 + // input[36]: defaulting = 0, visible = 0, type = 3, value = 1 + // input[41]: defaulting = 0, visible = 0, type = 32, value = "rotate" + // depth: value = 24 + // window: position = (0.1562,0.0426), size = 0.6641x0.8815 + // internal caching: 1 + // +main_Image_2_out_1, +main_Image_2_out_2, +main_Image_2_out_3 = + Image( + main_Image_2_in_1, + main_Collect_1_out_1, + main_Image_2_in_3, + main_Image_2_in_4, + main_Image_2_in_5, + main_Image_2_in_6, + main_Image_2_in_7, + main_Image_2_in_8, + main_Image_2_in_9, + main_Image_2_in_10, + main_Image_2_in_11, + main_Image_2_in_12, + main_Image_2_in_13, + main_Image_2_in_14, + main_Image_2_in_15, + main_Image_2_in_16, + main_Image_2_in_17, + main_Image_2_in_18, + main_Image_2_in_19, + main_Image_2_in_20, + main_Image_2_in_21, + main_Image_2_in_22, + main_Image_2_in_23, + main_Image_2_in_24, + main_Image_2_in_25, + main_Image_2_in_26, + main_Image_2_in_27, + main_Image_2_in_28, + main_Image_2_in_29, + main_Image_2_in_30, + main_Image_2_in_31, + main_Image_2_in_32, + main_Image_2_in_33, + main_Image_2_in_34, + main_Image_2_in_35, + main_Image_2_in_36, + main_Image_2_in_37, + main_Image_2_in_38, + main_Image_2_in_39, + main_Image_2_in_40, + main_Image_2_in_41, + main_Image_2_in_42, + main_Image_2_in_43, + main_Image_2_in_44, + main_Image_2_in_45, + main_Image_2_in_46, + main_Image_2_in_47, + main_Image_2_in_48, + main_Image_2_in_49 + ) [instance: 2, cache: 1]; +// network: end of macro body +CacheScene(main_Image_2_in_1, main_Image_2_out_1, main_Image_2_out_2); +} +main_AmbientLight_1_in_1 = NULL; +main_AmbientLight_1_out_1 = NULL; +main_Import_2_in_1 = "DX_fullmesh.dx"; +main_Import_2_in_2 = NULL; +main_Import_2_in_3 = "dx"; +main_Import_2_in_4 = NULL; +main_Import_2_in_5 = NULL; +main_Import_2_in_6 = NULL; +main_Import_2_out_1 = NULL; +main_AutoColor_1_in_2 = NULL; +main_AutoColor_1_in_3 = NULL; +main_AutoColor_1_in_4 = NULL; +main_AutoColor_1_in_5 = NULL; +main_AutoColor_1_in_6 = NULL; +main_AutoColor_1_in_7 = NULL; +main_AutoColor_1_in_8 = NULL; +main_AutoColor_1_in_9 = NULL; +main_AutoColor_1_in_10 = NULL; +main_AutoColor_1_out_1 = NULL; +main_ShowConnections_1_out_1 = NULL; +main_Collect_1_in_4 = NULL; +main_Collect_1_out_1 = NULL; +macro Image( + id, + object, + where, + useVector, + to, + from, + width, + resolution, + aspect, + up, + viewAngle, + perspective, + options, + buttonState = 1, + buttonUpApprox = "none", + buttonDownApprox = "none", + buttonUpDensity = 1, + buttonDownDensity = 1, + renderMode = 0, + defaultCamera, + reset, + backgroundColor, + throttle, + RECenable = 0, + RECfile, + RECformat, + RECresolution, + RECaspect, + AAenable = 0, + AAlabels, + AAticks, + AAcorners, + AAframe, + AAadjust, + AAcursor, + AAgrid, + AAcolors, + AAannotation, + AAlabelscale, + AAfont, + interactionMode, + title, + AAxTickLocs, + AAyTickLocs, + AAzTickLocs, + AAxTickLabels, + AAyTickLabels, + AAzTickLabels, + webOptions) -> ( + object, + camera, + where) +{ + ImageMessage( + id, + backgroundColor, + throttle, + RECenable, + RECfile, + RECformat, + RECresolution, + RECaspect, + AAenable, + AAlabels, + AAticks, + AAcorners, + AAframe, + AAadjust, + AAcursor, + AAgrid, + AAcolors, + AAannotation, + AAlabelscale, + AAfont, + AAxTickLocs, + AAyTickLocs, + AAzTickLocs, + AAxTickLabels, + AAyTickLabels, + AAzTickLabels, + interactionMode, + title, + renderMode, + buttonUpApprox, + buttonDownApprox, + buttonUpDensity, + buttonDownDensity) [instance: 1, cache: 1]; + autoCamera = + AutoCamera( + object, + "front", + object, + resolution, + aspect, + [0,1,0], + perspective, + viewAngle, + backgroundColor) [instance: 1, cache: 1]; + realCamera = + Camera( + to, + from, + width, + resolution, + aspect, + up, + perspective, + viewAngle, + backgroundColor) [instance: 1, cache: 1]; + coloredDefaultCamera = + UpdateCamera(defaultCamera, + background=backgroundColor) [instance: 1, cache: 1]; + nullDefaultCamera = + Inquire(defaultCamera, + "is null + 1") [instance: 1, cache: 1]; + resetCamera = + Switch( + nullDefaultCamera, + coloredDefaultCamera, + autoCamera) [instance: 1, cache: 1]; + resetNull = + Inquire( + reset, + "is null + 1") [instance: 2, cache: 1]; + reset = + Switch( + resetNull, + reset, + 0) [instance: 2, cache: 1]; + whichCamera = + Compute( + "($0 != 0 || $1 == 0) ? 1 : 2", + reset, + useVector) [instance: 1, cache: 1]; + camera = Switch( + whichCamera, + resetCamera, + realCamera) [instance: 3, cache: 1]; + AAobject = + AutoAxes( + object, + camera, + AAlabels, + AAticks, + AAcorners, + AAframe, + AAadjust, + AAcursor, + AAgrid, + AAcolors, + AAannotation, + AAlabelscale, + AAfont, + AAxTickLocs, + AAyTickLocs, + AAzTickLocs, + AAxTickLabels, + AAyTickLabels, + AAzTickLabels) [instance: 1, cache: 1]; + switchAAenable = Compute("$0+1", + AAenable) [instance: 2, cache: 1]; + object = Switch( + switchAAenable, + object, + AAobject) [instance:4, cache: 1]; + SWapproximation_options = + Switch( + buttonState, + buttonUpApprox, + buttonDownApprox) [instance: 5, cache: 1]; + SWdensity_options = + Switch( + buttonState, + buttonUpDensity, + buttonDownDensity) [instance: 6, cache: 1]; + HWapproximation_options = + Format( + "%s,%s", + buttonDownApprox, + buttonUpApprox) [instance: 1, cache: 1]; + HWdensity_options = + Format( + "%d,%d", + buttonDownDensity, + buttonUpDensity) [instance: 2, cache: 1]; + switchRenderMode = Compute( + "$0+1", + renderMode) [instance: 3, cache: 1]; + approximation_options = Switch( + switchRenderMode, + SWapproximation_options, + HWapproximation_options) [instance: 7, cache: 1]; + density_options = Switch( + switchRenderMode, + SWdensity_options, + HWdensity_options) [instance: 8, cache: 1]; + renderModeString = Switch( + switchRenderMode, + "software", + "hardware")[instance: 9, cache: 1]; + object_tag = Inquire( + object, + "object tag")[instance: 3, cache: 1]; + annoted_object = + Options( + object, + "send boxes", + 0, + "cache", + 1, + "object tag", + object_tag, + "ddcamera", + whichCamera, + "rendering approximation", + approximation_options, + "render every", + density_options, + "button state", + buttonState, + "rendering mode", + renderModeString) [instance: 1, cache: 1]; + RECresNull = + Inquire( + RECresolution, + "is null + 1") [instance: 4, cache: 1]; + ImageResolution = + Inquire( + camera, + "camera resolution") [instance: 5, cache: 1]; + RECresolution = + Switch( + RECresNull, + RECresolution, + ImageResolution) [instance: 10, cache: 1]; + RECaspectNull = + Inquire( + RECaspect, + "is null + 1") [instance: 6, cache: 1]; + ImageAspect = + Inquire( + camera, + "camera aspect") [instance: 7, cache: 1]; + RECaspect = + Switch( + RECaspectNull, + RECaspect, + ImageAspect) [instance: 11, cache: 1]; + switchRECenable = Compute( + "$0 == 0 ? 1 : (($2 == $3) && ($4 == $5)) ? ($1 == 1 ? 2 : 3) : 4", + RECenable, + switchRenderMode, + RECresolution, + ImageResolution, + RECaspect, + ImageAspect) [instance: 4, cache: 1]; + NoRECobject, RECNoRerenderObject, RECNoRerHW, RECRerenderObject = Route(switchRECenable, annoted_object); + Display( + NoRECobject, + camera, + where, + throttle) [instance: 1, cache: 1]; + image = + Render( + RECNoRerenderObject, + camera) [instance: 1, cache: 1]; + Display( + image, + NULL, + where, + throttle) [instance: 2, cache: 1]; + WriteImage( + image, + RECfile, + RECformat) [instance: 1, cache: 1]; + rec_where = Display( + RECNoRerHW, + camera, + where, + throttle) [instance: 1, cache: 0]; + rec_image = ReadImageWindow( + rec_where) [instance: 1, cache: 1]; + WriteImage( + rec_image, + RECfile, + RECformat) [instance: 1, cache: 1]; + RECupdateCamera = + UpdateCamera( + camera, + resolution=RECresolution, + aspect=RECaspect) [instance: 2, cache: 1]; + Display( + RECRerenderObject, + camera, + where, + throttle) [instance: 1, cache: 1]; + RECRerenderObject = + ScaleScreen( + RECRerenderObject, + NULL, + RECresolution, + camera) [instance: 1, cache: 1]; + image = + Render( + RECRerenderObject, + RECupdateCamera) [instance: 2, cache: 1]; + WriteImage( + image, + RECfile, + RECformat) [instance: 2, cache: 1]; +} +main_Image_2_in_1 = "Image_2"; +main_Image_2_in_3 = "X24,,"; +main_Image_2_in_4 = 1; +main_Image_2_in_5 = [0.025 0.025 0.0425]; +main_Image_2_in_6 = [0.0943607 -0.252846 0.0756326]; +main_Image_2_in_7 = 0.154491; +main_Image_2_in_8 = 1261; +main_Image_2_in_9 = 0.722839; +main_Image_2_in_10 = [0.0387653 0.127856 0.991035]; +main_Image_2_in_11 = NULL; +main_Image_2_in_12 = 0; +main_Image_2_in_13 = NULL; +main_Image_2_in_14 = 1; +main_Image_2_in_15 = "none"; +main_Image_2_in_16 = "none"; +main_Image_2_in_17 = 1; +main_Image_2_in_18 = 1; +main_Image_2_in_19 = 1; +main_Image_2_in_20 = NULL; +main_Image_2_in_21 = NULL; +main_Image_2_in_22 = "black"; +main_Image_2_in_23 = NULL; +main_Image_2_in_25 = "/home/latychev/boo"; +main_Image_2_in_26 = "eps gray dpi=86 orient=portrait"; +main_Image_2_in_27 = NULL; +main_Image_2_in_28 = NULL; +main_Image_2_in_29 = 0; +main_Image_2_in_30 = NULL; +main_Image_2_in_31 = NULL; +main_Image_2_in_32 = NULL; +main_Image_2_in_33 = 1; +main_Image_2_in_34 = NULL; +main_Image_2_in_35 = NULL; +main_Image_2_in_36 = 1; +main_Image_2_in_37 = NULL; +main_Image_2_in_38 = NULL; +main_Image_2_in_39 = NULL; +main_Image_2_in_40 = NULL; +main_Image_2_in_41 = "rotate"; +main_Image_2_in_42 = NULL; +main_Image_2_in_43 = NULL; +main_Image_2_in_44 = NULL; +main_Image_2_in_45 = NULL; +main_Image_2_in_46 = NULL; +main_Image_2_in_47 = NULL; +main_Image_2_in_48 = NULL; +main_Image_2_in_49 = NULL; +Executive("product version 4 4 4"); +$sync +main();