[cig-commits] r15941 - seismo/2D/SPECFEM2D/trunk

[cmorency at geodynamics.org]

Author: cmorency
Date: 2009-11-09 12:51:21 -0800 (Mon, 09 Nov 2009)
New Revision: 15941

Modified:
   seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
Log:
Modified detection of free surfaces for acoustic case and regular mesh. Only the top surface was detected. Now all the 4 surfaces are set to free unless absorbing condition is prescribed.  


Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90	2009-11-09 17:39:55 UTC (rev 15940)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90	2009-11-09 20:51:21 UTC (rev 15941)
@@ -1091,9 +1091,8 @@
      ! count the number of acoustic free-surface elements
      nelem_acoustic_surface = 0
 
-! if the top surface is absorbing, it cannot be free at the same time
+! if the surface is absorbing, it cannot be free at the same time
      if(.not. abstop) then
-
      j = nzread
      do i = 1,nxread
         imaterial_number = num_material((j-1)*nxread+i)
@@ -1101,10 +1100,41 @@
            nelem_acoustic_surface = nelem_acoustic_surface + 1
         endif
      enddo
+     endif
+     if(.not. absbottom) then
+     j = 1
+     do i = 1,nxread
+        imaterial_number = num_material((j-1)*nxread+i)
+        if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
+           nelem_acoustic_surface = nelem_acoustic_surface + 1
+        endif
+     enddo
+     endif
+     if(.not. absleft) then
+     i = 1
+     do j = 1,nzread
+        imaterial_number = num_material((j-1)*nxread+i)
+        if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
+           nelem_acoustic_surface = nelem_acoustic_surface + 1
+        endif
+     enddo
+     endif
+     if(.not. absright) then
+     i = nxread
+     do j = 1,nzread
+        imaterial_number = num_material((j-1)*nxread+i)
+        if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
+           nelem_acoustic_surface = nelem_acoustic_surface + 1
+        endif
+     enddo
+     endif
 
+
      allocate(acoustic_surface(4,nelem_acoustic_surface))
 
      nelem_acoustic_surface = 0
+
+     if(.not. abstop) then
      j = nzread
      do i = 1,nxread
         imaterial_number = num_material((j-1)*nxread+i)
@@ -1116,8 +1146,46 @@
            acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
         endif
      enddo
-
      endif
+     if(.not. absbottom) then
+     j = 1
+     do i = 1,nxread
+        imaterial_number = num_material((j-1)*nxread+i)
+        if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
+           nelem_acoustic_surface = nelem_acoustic_surface + 1
+           acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
+           acoustic_surface(2,nelem_acoustic_surface) = 2
+           acoustic_surface(3,nelem_acoustic_surface) = elmnts(0+ngnod*((j-1)*nxread+i-1))
+           acoustic_surface(4,nelem_acoustic_surface) = elmnts(1+ngnod*((j-1)*nxread+i-1))
+        endif
+     enddo
+     endif
+     if(.not. absleft) then
+     i = 1
+     do j = 1,nzread
+        imaterial_number = num_material((j-1)*nxread+i)
+        if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
+           nelem_acoustic_surface = nelem_acoustic_surface + 1
+           acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
+           acoustic_surface(2,nelem_acoustic_surface) = 2
+           acoustic_surface(3,nelem_acoustic_surface) = elmnts(0+ngnod*((j-1)*nxread+i-1))
+           acoustic_surface(4,nelem_acoustic_surface) = elmnts(3+ngnod*((j-1)*nxread+i-1))
+        endif
+     enddo
+     endif
+     if(.not. absright) then
+     i = nxread
+     do j = 1,nzread
+        imaterial_number = num_material((j-1)*nxread+i)
+        if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
+           nelem_acoustic_surface = nelem_acoustic_surface + 1
+           acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
+           acoustic_surface(2,nelem_acoustic_surface) = 2
+           acoustic_surface(3,nelem_acoustic_surface) = elmnts(1+ngnod*((j-1)*nxread+i-1))
+           acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
+        endif
+     enddo
+     endif
 
      !
      !--- definition of absorbing boundaries