[cig-commits] [commit] devel: Modify the file list_for_Zhinan_for_SPECFEM2D_PMLs.txt, since the plane incident bug have already been fixed by Dimitri. See the subroutine rotate_mesh_for_plane_wave in part_unstruct.F90. That subroutine is adapted from two of Dimitri's old codes, which also added again in UTILS/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries. (f7f8aed)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Jan 8 01:22:44 PST 2015
Repository : https://github.com/geodynamics/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/de54f292fc495b3e612e883ea9c3ad063282a6e9...f7f8aedbad5f120935b18e9c1d117c81aa19aa06
>---------------------------------------------------------------
commit f7f8aedbad5f120935b18e9c1d117c81aa19aa06
Author: Xie Zhinan <xiezhinan1984 at gmail.com>
Date: Thu Jan 8 10:03:04 2015 +0800
Modify the file list_for_Zhinan_for_SPECFEM2D_PMLs.txt,
since the plane incident bug have already been fixed by Dimitri.
See the subroutine rotate_mesh_for_plane_wave in part_unstruct.F90.
That subroutine is adapted from two of Dimitri's old codes,
which also added again in UTILS/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries.
>---------------------------------------------------------------
f7f8aedbad5f120935b18e9c1d117c81aa19aa06
..._of_external_mesh_for_plane_wave_boundaries.f90 | 154 ++++++++
.../rotate_mesh_for_absorbing_conditions.f90 | 426 +++++++++++++++++++++
list_for_Zhinan_for_SPECFEM2D_PMLs.txt | 60 ---
3 files changed, 580 insertions(+), 60 deletions(-)
diff --git a/UTILS/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries.f90 b/UTILS/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries.f90
new file mode 100644
index 0000000..b32723b
--- /dev/null
+++ b/UTILS/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries.f90
@@ -0,0 +1,154 @@
+ program djjkfdkjdf
+
+implicit none
+
+integer, parameter :: ngnod = 4 !! DK DK for now, will switch to 9 later
+
+double precision, parameter :: x_left_edge = -1.d0 !! DK DK for now
+double precision, parameter :: x_right_edge = +1.d0 !! DK DK for now
+
+double precision, parameter :: tol = (x_right_edge - x_left_edge) / 1.d6
+
+double precision, dimension(:), allocatable :: x,z
+integer, dimension(:,:), allocatable :: ibool
+
+integer :: npoin,nspec,ipoin,ispec,counter
+
+! ========================== first process the left edge ==========================
+
+! header
+open(unit=12,file='Nodes_SqrCircMeta.msh',status='old')
+
+read(12,*) npoin
+allocate(x(npoin))
+allocate(z(npoin))
+
+do ipoin = 1,npoin
+ read(12,*) x(ipoin),z(ipoin)
+enddo
+close(12)
+
+! header
+open(unit=12,file='Mesh_SqrCircMeta.msh',status='old')
+
+read(12,*) nspec
+allocate(ibool(ngnod,nspec))
+
+do ispec = 1,nspec
+ read(12,*) ibool(1,ispec),ibool(2,ispec),ibool(3,ispec),ibool(4,ispec)
+enddo
+
+close(12)
+
+open(unit=12,file='Mesh_SqrCircMeta.msh',status='unknown')
+
+write(12,*) nspec
+
+! find left edge
+counter = 0
+do ispec = 1,nspec
+
+ if(abs(x(ibool(1,ispec)) - x_left_edge) < tol .and. abs(x(ibool(2,ispec)) - x_left_edge) < tol) then
+ print *,'elem ',ispec,' is on left edge with edge 1'
+ counter = counter + 1
+ write(12,*) ibool(2,ispec),ibool(3,ispec),ibool(4,ispec),ibool(1,ispec)
+
+ else if(abs(x(ibool(2,ispec)) - x_left_edge) < tol .and. abs(x(ibool(3,ispec)) - x_left_edge) < tol) then
+ print *,'elem ',ispec,' is on left edge with edge 2'
+ counter = counter + 1
+ write(12,*) ibool(3,ispec),ibool(4,ispec),ibool(1,ispec),ibool(2,ispec)
+
+ else if(abs(x(ibool(3,ispec)) - x_left_edge) < tol .and. abs(x(ibool(4,ispec)) - x_left_edge) < tol) then
+ print *,'elem ',ispec,' is on left edge with edge 3'
+ counter = counter + 1
+ write(12,*) ibool(4,ispec),ibool(1,ispec),ibool(2,ispec),ibool(3,ispec)
+
+ else if(abs(x(ibool(4,ispec)) - x_left_edge) < tol .and. abs(x(ibool(1,ispec)) - x_left_edge) < tol) then
+ print *,'elem ',ispec,' is on left edge with edge 4'
+ counter = counter + 1
+ write(12,*) ibool(1,ispec),ibool(2,ispec),ibool(3,ispec),ibool(4,ispec)
+
+ else
+ write(12,*) ibool(1,ispec),ibool(2,ispec),ibool(3,ispec),ibool(4,ispec)
+ endif
+
+enddo
+
+print *
+print *,'found ',counter,' elements on left edge'
+print *
+
+close(12)
+
+deallocate(x,z,ibool)
+
+! ========================== then process the right edge ==========================
+
+! header
+open(unit=12,file='Nodes_SqrCircMeta.msh',status='old')
+
+read(12,*) npoin
+allocate(x(npoin))
+allocate(z(npoin))
+
+do ipoin = 1,npoin
+ read(12,*) x(ipoin),z(ipoin)
+enddo
+close(12)
+
+! header
+open(unit=12,file='Mesh_SqrCircMeta.msh',status='old')
+
+read(12,*) nspec
+allocate(ibool(ngnod,nspec))
+
+do ispec = 1,nspec
+ read(12,*) ibool(1,ispec),ibool(2,ispec),ibool(3,ispec),ibool(4,ispec)
+enddo
+
+close(12)
+
+open(unit=12,file='Mesh_SqrCircMeta.msh',status='unknown')
+
+write(12,*) nspec
+
+! find right edge
+counter = 0
+do ispec = 1,nspec
+
+ if(abs(x(ibool(1,ispec)) - x_right_edge) < tol .and. abs(x(ibool(2,ispec)) - x_right_edge) < tol) then
+ print *,'elem ',ispec,' is on right edge with edge 1'
+ counter = counter + 1
+ write(12,*) ibool(4,ispec),ibool(1,ispec),ibool(2,ispec),ibool(3,ispec)
+
+ else if(abs(x(ibool(2,ispec)) - x_right_edge) < tol .and. abs(x(ibool(3,ispec)) - x_right_edge) < tol) then
+ print *,'elem ',ispec,' is on right edge with edge 2'
+ counter = counter + 1
+ write(12,*) ibool(1,ispec),ibool(2,ispec),ibool(3,ispec),ibool(4,ispec)
+
+ else if(abs(x(ibool(3,ispec)) - x_right_edge) < tol .and. abs(x(ibool(4,ispec)) - x_right_edge) < tol) then
+ print *,'elem ',ispec,' is on right edge with edge 3'
+ counter = counter + 1
+ write(12,*) ibool(2,ispec),ibool(3,ispec),ibool(4,ispec),ibool(1,ispec)
+
+ else if(abs(x(ibool(4,ispec)) - x_right_edge) < tol .and. abs(x(ibool(1,ispec)) - x_right_edge) < tol) then
+ print *,'elem ',ispec,' is on right edge with edge 4'
+ counter = counter + 1
+ write(12,*) ibool(3,ispec),ibool(4,ispec),ibool(1,ispec),ibool(2,ispec)
+
+ else
+ write(12,*) ibool(1,ispec),ibool(2,ispec),ibool(3,ispec),ibool(4,ispec)
+ endif
+
+enddo
+
+print *
+print *,'found ',counter,' elements on right edge'
+print *
+
+close(12)
+
+deallocate(x,z,ibool)
+
+end
+
diff --git a/UTILS/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries/rotate_mesh_for_absorbing_conditions.f90 b/UTILS/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries/rotate_mesh_for_absorbing_conditions.f90
new file mode 100644
index 0000000..b999d66
--- /dev/null
+++ b/UTILS/program_rotate_edges_of_external_mesh_for_plane_wave_boundaries/rotate_mesh_for_absorbing_conditions.f90
@@ -0,0 +1,426 @@
+
+ program rotate_mesh
+
+ implicit none
+
+!! DK DK Dec 2011: rotate mesh elements to make sure topological absorbing surface is aligned with physical absorbing surface
+
+integer, parameter :: ngnod = 9
+
+integer :: nspec, nglob, nabs
+
+integer :: i,ispec,iglob,idummy,i1,i2,inode,iswap
+
+logical :: found_this_point
+
+integer, dimension(:,:), allocatable :: ibool,ibool_rotated
+
+integer, dimension(:), allocatable :: ielem,iglob1,iglob2
+
+open(unit=12,file='Mesh_ht_test1',status='old')
+read(12,*) nspec
+allocate(ibool(ngnod,nspec))
+allocate(ibool_rotated(ngnod,nspec))
+do ispec = 1,nspec
+ read(12,*) ibool(1,ispec), ibool(2,ispec), ibool(3,ispec), ibool(4,ispec), &
+ ibool(5,ispec), ibool(6,ispec), ibool(7,ispec), ibool(8,ispec), ibool(9,ispec)
+enddo
+close(12)
+
+!!!!!!!!!!! rotate for TOP surface !!!!!!!!!!!!!
+
+! first make a copy of ibool (for all elements that will be unchanged)
+ ibool_rotated(:,:) = ibool(:,:)
+
+open(unit=12,file='Surf_top_ht_test1',status='old')
+read(12,*) nabs
+allocate(ielem(nabs))
+allocate(iglob1(nabs))
+allocate(iglob2(nabs))
+do i = 1,nabs
+ read(12,*) ielem(i),idummy,iglob1(i),iglob2(i)
+enddo
+close(12)
+
+! look for rotations needed
+do i = 1,nabs
+ ispec = ielem(i)
+
+ found_this_point = .false.
+ do inode = 1,4
+ if(ibool(inode,ispec) == iglob1(i)) then
+ i1 = inode
+ found_this_point = .true.
+ exit
+ endif
+ enddo
+ if(.not. found_this_point) stop 'point not found'
+
+ found_this_point = .false.
+ do inode = 1,4
+ if(ibool(inode,ispec) == iglob2(i)) then
+ i2 = inode
+ found_this_point = .true.
+ exit
+ endif
+ enddo
+ if(.not. found_this_point) stop 'point not found'
+
+! swap points if needed for clarity, to avoid testing both cases each time below
+ if(i1 > i2) then
+ iswap = i1
+ i1 = i2
+ i2 = iswap
+ endif
+
+! test orientation
+ if(i1 == 3 .and. i2 == 4) then
+ print *,'orientation of element ',i,' is already good'
+
+ else if (i1 == 1 .and. i2 == 4) then ! for this one, remember that we have swapped, thus 4 1 is 1 4
+ print *,'element ',i,' must be rotated 3 times'
+ ibool_rotated(4,ispec) = ibool(1,ispec)
+ ibool_rotated(1,ispec) = ibool(2,ispec)
+ ibool_rotated(2,ispec) = ibool(3,ispec)
+ ibool_rotated(3,ispec) = ibool(4,ispec)
+ ibool_rotated(8,ispec) = ibool(5,ispec)
+ ibool_rotated(5,ispec) = ibool(6,ispec)
+ ibool_rotated(6,ispec) = ibool(7,ispec)
+ ibool_rotated(7,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else if (i1 == 1 .and. i2 == 2) then
+ print *,'element ',i,ispec,' must be rotated 2 times top'
+ ibool_rotated(3,ispec) = ibool(1,ispec)
+ ibool_rotated(4,ispec) = ibool(2,ispec)
+ ibool_rotated(1,ispec) = ibool(3,ispec)
+ ibool_rotated(2,ispec) = ibool(4,ispec)
+ ibool_rotated(7,ispec) = ibool(5,ispec)
+ ibool_rotated(8,ispec) = ibool(6,ispec)
+ ibool_rotated(5,ispec) = ibool(7,ispec)
+ ibool_rotated(6,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else if (i1 == 2 .and. i2 == 3) then
+ print *,'element ',i,' must be rotated 1 time'
+ ibool_rotated(2,ispec) = ibool(1,ispec)
+ ibool_rotated(3,ispec) = ibool(2,ispec)
+ ibool_rotated(4,ispec) = ibool(3,ispec)
+ ibool_rotated(1,ispec) = ibool(4,ispec)
+ ibool_rotated(6,ispec) = ibool(5,ispec)
+ ibool_rotated(7,ispec) = ibool(6,ispec)
+ ibool_rotated(8,ispec) = ibool(7,ispec)
+ ibool_rotated(5,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else
+ stop 'problem in an element'
+ endif
+
+enddo
+
+deallocate(ielem)
+deallocate(iglob1)
+deallocate(iglob2)
+
+!!!!!!!!!!! rotate for BOTTOM surface !!!!!!!!!!!!!
+
+! copy the changes above back to ibool, which is used again below
+ ibool(:,:) = ibool_rotated(:,:)
+
+open(unit=12,file='Surf_bottom_ht_test1',status='old')
+read(12,*) nabs
+allocate(ielem(nabs))
+allocate(iglob1(nabs))
+allocate(iglob2(nabs))
+do i = 1,nabs
+ read(12,*) ielem(i),idummy,iglob1(i),iglob2(i)
+enddo
+close(12)
+
+! look for rotations needed
+do i = 1,nabs
+ ispec = ielem(i)
+
+ found_this_point = .false.
+ do inode = 1,4
+ if(ibool(inode,ispec) == iglob1(i)) then
+ i1 = inode
+ found_this_point = .true.
+ exit
+ endif
+ enddo
+ if(.not. found_this_point) stop 'point not found'
+
+ found_this_point = .false.
+ do inode = 1,4
+ if(ibool(inode,ispec) == iglob2(i)) then
+ i2 = inode
+ found_this_point = .true.
+ exit
+ endif
+ enddo
+ if(.not. found_this_point) stop 'point not found'
+
+! swap points if needed for clarity, to avoid testing both cases each time below
+ if(i1 > i2) then
+ iswap = i1
+ i1 = i2
+ i2 = iswap
+ endif
+
+! test orientation
+ if(i1 == 1 .and. i2 == 2) then
+ print *,'orientation of element ',i,' is already good'
+
+ else if (i1 == 2 .and. i2 == 3) then
+ print *,'element ',i,' must be rotated 3 times'
+ ibool_rotated(4,ispec) = ibool(1,ispec)
+ ibool_rotated(1,ispec) = ibool(2,ispec)
+ ibool_rotated(2,ispec) = ibool(3,ispec)
+ ibool_rotated(3,ispec) = ibool(4,ispec)
+ ibool_rotated(8,ispec) = ibool(5,ispec)
+ ibool_rotated(5,ispec) = ibool(6,ispec)
+ ibool_rotated(6,ispec) = ibool(7,ispec)
+ ibool_rotated(7,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else if (i1 == 3 .and. i2 == 4) then
+ print *,'element ',i,ispec,' must be rotated 2 times bottom'
+ ibool_rotated(3,ispec) = ibool(1,ispec)
+ ibool_rotated(4,ispec) = ibool(2,ispec)
+ ibool_rotated(1,ispec) = ibool(3,ispec)
+ ibool_rotated(2,ispec) = ibool(4,ispec)
+ ibool_rotated(7,ispec) = ibool(5,ispec)
+ ibool_rotated(8,ispec) = ibool(6,ispec)
+ ibool_rotated(5,ispec) = ibool(7,ispec)
+ ibool_rotated(6,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else if (i1 == 1 .and. i2 == 4) then ! for this one, remember that we have swapped, thus 4 1 is 1 4
+ print *,'element ',i,' must be rotated 1 time'
+ ibool_rotated(2,ispec) = ibool(1,ispec)
+ ibool_rotated(3,ispec) = ibool(2,ispec)
+ ibool_rotated(4,ispec) = ibool(3,ispec)
+ ibool_rotated(1,ispec) = ibool(4,ispec)
+ ibool_rotated(6,ispec) = ibool(5,ispec)
+ ibool_rotated(7,ispec) = ibool(6,ispec)
+ ibool_rotated(8,ispec) = ibool(7,ispec)
+ ibool_rotated(5,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else
+ stop 'problem in an element'
+ endif
+
+enddo
+
+deallocate(ielem)
+deallocate(iglob1)
+deallocate(iglob2)
+
+!!!!!!!!!!! rotate for LEFT surface !!!!!!!!!!!!!
+
+! copy the changes above back to ibool, which is used again below
+ ibool(:,:) = ibool_rotated(:,:)
+
+open(unit=12,file='Surf_left_ht_test1',status='old')
+read(12,*) nabs
+allocate(ielem(nabs))
+allocate(iglob1(nabs))
+allocate(iglob2(nabs))
+do i = 1,nabs
+ read(12,*) ielem(i),idummy,iglob1(i),iglob2(i)
+enddo
+close(12)
+
+! look for rotations needed
+do i = 1,nabs
+ ispec = ielem(i)
+
+ found_this_point = .false.
+ do inode = 1,4
+ if(ibool(inode,ispec) == iglob1(i)) then
+ i1 = inode
+ found_this_point = .true.
+ exit
+ endif
+ enddo
+ if(.not. found_this_point) stop 'point not found'
+
+ found_this_point = .false.
+ do inode = 1,4
+ if(ibool(inode,ispec) == iglob2(i)) then
+ i2 = inode
+ found_this_point = .true.
+ exit
+ endif
+ enddo
+ if(.not. found_this_point) stop 'point not found'
+
+! swap points if needed for clarity, to avoid testing both cases each time below
+ if(i1 > i2) then
+ iswap = i1
+ i1 = i2
+ i2 = iswap
+ endif
+
+! test orientation
+ if(i1 == 1 .and. i2 == 4) then ! for this one, remember that we have swapped, thus 4 1 is 1 4
+ print *,'orientation of element ',i,' is already good'
+
+ else if (i1 == 1 .and. i2 == 2) then
+ print *,'element ',i,' must be rotated 3 times'
+ ibool_rotated(4,ispec) = ibool(1,ispec)
+ ibool_rotated(1,ispec) = ibool(2,ispec)
+ ibool_rotated(2,ispec) = ibool(3,ispec)
+ ibool_rotated(3,ispec) = ibool(4,ispec)
+ ibool_rotated(8,ispec) = ibool(5,ispec)
+ ibool_rotated(5,ispec) = ibool(6,ispec)
+ ibool_rotated(6,ispec) = ibool(7,ispec)
+ ibool_rotated(7,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else if (i1 == 2 .and. i2 == 3) then
+ print *,'element ',i,ispec,' must be rotated 2 times left'
+ ibool_rotated(3,ispec) = ibool(1,ispec)
+ ibool_rotated(4,ispec) = ibool(2,ispec)
+ ibool_rotated(1,ispec) = ibool(3,ispec)
+ ibool_rotated(2,ispec) = ibool(4,ispec)
+ ibool_rotated(7,ispec) = ibool(5,ispec)
+ ibool_rotated(8,ispec) = ibool(6,ispec)
+ ibool_rotated(5,ispec) = ibool(7,ispec)
+ ibool_rotated(6,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else if (i1 == 3 .and. i2 == 4) then
+ print *,'element ',i,' must be rotated 1 time'
+ ibool_rotated(2,ispec) = ibool(1,ispec)
+ ibool_rotated(3,ispec) = ibool(2,ispec)
+ ibool_rotated(4,ispec) = ibool(3,ispec)
+ ibool_rotated(1,ispec) = ibool(4,ispec)
+ ibool_rotated(6,ispec) = ibool(5,ispec)
+ ibool_rotated(7,ispec) = ibool(6,ispec)
+ ibool_rotated(8,ispec) = ibool(7,ispec)
+ ibool_rotated(5,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else
+ stop 'problem in an element'
+ endif
+
+enddo
+
+deallocate(ielem)
+deallocate(iglob1)
+deallocate(iglob2)
+
+!!!!!!!!!!! rotate for RIGHT surface !!!!!!!!!!!!!
+
+! copy the changes above back to ibool, which is used again below
+ ibool(:,:) = ibool_rotated(:,:)
+
+open(unit=12,file='Surf_right_ht_test1',status='old')
+read(12,*) nabs
+allocate(ielem(nabs))
+allocate(iglob1(nabs))
+allocate(iglob2(nabs))
+do i = 1,nabs
+ read(12,*) ielem(i),idummy,iglob1(i),iglob2(i)
+enddo
+close(12)
+
+! look for rotations needed
+do i = 1,nabs
+ ispec = ielem(i)
+
+ found_this_point = .false.
+ do inode = 1,4
+ if(ibool(inode,ispec) == iglob1(i)) then
+ i1 = inode
+ found_this_point = .true.
+ exit
+ endif
+ enddo
+ if(.not. found_this_point) stop 'point not found'
+
+ found_this_point = .false.
+ do inode = 1,4
+ if(ibool(inode,ispec) == iglob2(i)) then
+ i2 = inode
+ found_this_point = .true.
+ exit
+ endif
+ enddo
+ if(.not. found_this_point) stop 'point not found'
+
+! swap points if needed for clarity, to avoid testing both cases each time below
+ if(i1 > i2) then
+ iswap = i1
+ i1 = i2
+ i2 = iswap
+ endif
+
+! test orientation
+ if(i1 == 2 .and. i2 == 3) then
+ print *,'orientation of element ',i,' is already good'
+
+ else if (i1 == 3 .and. i2 == 4) then
+ print *,'element ',i,' must be rotated 3 times'
+ ibool_rotated(4,ispec) = ibool(1,ispec)
+ ibool_rotated(1,ispec) = ibool(2,ispec)
+ ibool_rotated(2,ispec) = ibool(3,ispec)
+ ibool_rotated(3,ispec) = ibool(4,ispec)
+ ibool_rotated(8,ispec) = ibool(5,ispec)
+ ibool_rotated(5,ispec) = ibool(6,ispec)
+ ibool_rotated(6,ispec) = ibool(7,ispec)
+ ibool_rotated(7,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else if (i1 == 1 .and. i2 == 4) then ! for this one, remember that we have swapped, thus 4 1 is 1 4
+ print *,'element ',i,ispec,' must be rotated 2 times right'
+ ibool_rotated(3,ispec) = ibool(1,ispec)
+ ibool_rotated(4,ispec) = ibool(2,ispec)
+ ibool_rotated(1,ispec) = ibool(3,ispec)
+ ibool_rotated(2,ispec) = ibool(4,ispec)
+ ibool_rotated(7,ispec) = ibool(5,ispec)
+ ibool_rotated(8,ispec) = ibool(6,ispec)
+ ibool_rotated(5,ispec) = ibool(7,ispec)
+ ibool_rotated(6,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else if (i1 == 1 .and. i2 == 2) then
+ print *,'element ',i,' must be rotated 1 time'
+ ibool_rotated(2,ispec) = ibool(1,ispec)
+ ibool_rotated(3,ispec) = ibool(2,ispec)
+ ibool_rotated(4,ispec) = ibool(3,ispec)
+ ibool_rotated(1,ispec) = ibool(4,ispec)
+ ibool_rotated(6,ispec) = ibool(5,ispec)
+ ibool_rotated(7,ispec) = ibool(6,ispec)
+ ibool_rotated(8,ispec) = ibool(7,ispec)
+ ibool_rotated(5,ispec) = ibool(8,ispec)
+! 9th point is at the element center and thus never changes when we rotate an element
+
+ else
+ stop 'problem in an element'
+ endif
+
+enddo
+
+deallocate(ielem)
+deallocate(iglob1)
+deallocate(iglob2)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! save the rotated mesh file
+open(unit=12,file='Mesh_ht_test1_rotated',status='unknown')
+write(12,*) nspec
+do ispec = 1,nspec
+ write(12,*) ibool_rotated(1,ispec), ibool_rotated(2,ispec), ibool_rotated(3,ispec), ibool_rotated(4,ispec), &
+ ibool_rotated(5,ispec), ibool_rotated(6,ispec), ibool_rotated(7,ispec), ibool_rotated(8,ispec), ibool_rotated(9,ispec)
+enddo
+close(12)
+
+end
+
diff --git a/list_for_Zhinan_for_SPECFEM2D_PMLs.txt b/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
index 148482a..098d33f 100644
--- a/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
+++ b/list_for_Zhinan_for_SPECFEM2D_PMLs.txt
@@ -14,63 +14,3 @@ currently, we still use potential as the adjoint source, we need to change to pr
put a flags in constants.h.in that decide do PML_parameter_adjustment implemented in pml_init.F90 or not
-------------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------------
-
-5/ plane wave boundary conditions maybe have a (small) bug in SPECFEM2D, see https://github.com/geodynamics/specfem2d/issues/83
-
-On 07/03/2014 09:02, è°¢å¿—å— wrote:
-> Dear Dimitri,
->
-> I think this bug we have not fixed yet.
-> I will fix this during the commit of kernel code.
->
-> I am really sorry for the delay.
->
-> Thank you so much.
-> Best regards,
-> Zhinan
->
->
-> 2014-03-05 2:45 GMT+08:00 Dimitri Komatitsch:
->
-> Dear Zhinan,
->
-> Do you remember if this bug (plane waves that did not work fine any
-> more on the edges of the grid in the 2D code) has been fixed or not?
->
-> (I think so, but I am not 100% sure; if so, please let me know, and
-> I will then remove it from the list of things to do).
->
-> Thank you very much,
-> Best regards,
->
-> Dimitri.
->
-> On 18/09/2013 16:50, Dimitri Komatitsch wrote:
->
-> Dear Diego,
->
-> Zhinan has told me that he will fix this next week.
-> Sorry about the delay, we were finishing another project. He
-> will send
-> you an email when the problem is fixed and will cc us.
->
-> Thanks,
-> Cheers,
->
-> Dimitri.
->
-> On 05/06/2013 12:32 AM, Dimitri Komatitsch wrote:
->>
-> Dear Diego,
->
-> Florian Cachoux (a student of Raphael Garcia and me) has
-> found the same
-> bug. Zhinan is currently investigating. The bug will be
-> fixed soon
-> (sometime this month). We will let you know. I cc everybody.
->
-> Thanks,
-> Cheers,
->
-> Dimitri.
-
More information about the CIG-COMMITS
mailing list