[cig-commits] r19324 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES EXAMPLES/Rayleigh_wave_no_crack EXAMPLES/Rayleigh_wave_with_crack src/meshfem2D src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Dec 28 16:40:50 PST 2011
Author: dkomati1
Date: 2011-12-28 16:40:50 -0800 (Wed, 28 Dec 2011)
New Revision: 19324
Added:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack/
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/Par_file_Rayleigh_2D
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/README
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/SOURCE_Rayleigh_2D
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/interfaces_Rayleigh_flat.dat
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/process.sh
Removed:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/
Modified:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack/README
seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_regions.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
added an example with an incident Rayleigh wave propagating towards a crack; modified the source code in order to "manually" add a crack (i.e., some parameters are hardwired and the crack detection and handling code I added is not very flexible yet)
Copied: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack (from rev 19323, seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave)
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack/README
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/README 2011-12-28 14:36:37 UTC (rev 19323)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_no_crack/README 2011-12-29 00:40:50 UTC (rev 19324)
@@ -14,7 +14,7 @@
2. compile:
make all
-3. cd EXAMPLES/Rayleigh_wave
+3. cd EXAMPLES/Rayleigh_wave_no_crack
4. execute script to run mesher and solver for the PSV case:
./process.sh
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/Par_file_Rayleigh_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/Par_file_Rayleigh_2D (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/Par_file_Rayleigh_2D 2011-12-29 00:40:50 UTC (rev 19324)
@@ -0,0 +1,139 @@
+# title of job
+title = Rayleigh wave in a homogeneous medium with crack
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
+NOISE_TOMOGRAPHY = 0 # 0 = earthquake simulation, 1/2/3 = noise simulation
+SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
+
+# parameters concerning partitioning
+nproc = 1 # number of processes
+partitioning_method = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE = .false. # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod = 4 # number of control nodes per element (4 or 9)
+initialfield = .true. # use a plane wave as source or not
+add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
+assign_external_model = .false. # define external earth model or not
+READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
+TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
+Q0 = 1 # quality factor for viscous attenuation
+freq0 = 10 # frequency for viscous attenuation
+p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt = 5000 # total number of time steps
+deltat = 3.90625e-3 # duration of a time step
+USER_T0 = 0.0d0 # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS = 2 # number of standard linear solids for attenuation
+f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver line parameters for seismograms
+seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS = .true. # creates a STATION file in ./DATA
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT = .false. # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver line
+nrec = 5 # number of receivers
+xdeb = 14.5 # first receiver x in meters
+zdeb = 9. # first receiver z in meters
+xfin = 17.5 # last receiver x in meters (ignored if onlyone receiver)
+zfin = 9. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO = 200 # display frequency in time steps
+output_postscript_snapshot = .true. # output Postscript snapshot of the results
+output_color_image = .true. # output color image of the results
+imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps = 1. # minimum amplitude in % for snapshots
+meshvect = .true. # display mesh on vector plots or not
+modelvect = .false. # display velocity model on vector plots
+boundvect = .true. # display boundary conditions on plots
+interpol = .true. # interpolation of the display or not
+pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript = 1 # subsampling of background velocity model in PostScript snapshots
+factor_subsample_image = 1 # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR = 0.30d0 # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true. # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
+US_LETTER = .false. # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false. # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot = .false. # generate a GNUPLOT file for the grid
+output_grid = .false. # save the grid in a text file or not
+output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels = 3 # nb of different models
+# define models as
+# I: (model_number 1 rho Vp Vs 0 0 QKappa Qmu 0 0 0 0 0 0) or
+# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 1.d0 2.d0 1.d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
+## DK DK create three identical media (thus it is still a homogeneous medium)
+## DK DK just to have a way of detecting the crack in the mesh and to indicate with colors where the crack is located in the PostScript displays
+2 1 1.d0 2.d0 1.d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
+3 1 1.d0 2.d0 1.d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh = .false.
+
+# absorbing boundary active or not
+absorbing_conditions = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file = ./DATA/mesh/canyon_mesh_file # file containing the mesh
+nodes_coords_file = ./DATA/mesh/canyon_nodes_coords_file # file containing the nodes coordinates
+materials_file = ./DATA/mesh/canyon_materials_file # file containing the material number for each element
+free_surface_file = ./DATA/mesh/canyon_free_surface_file # file containing the free surface
+absorbing_surface_file = ./DATA/mesh/canyon_absorbing_surface_file # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile = interfaces_rayleigh_flat.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 19.d0 # abscissa of right side of the model
+nx = 60 # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom = .true.
+absorbright = .true.
+absorbtop = .false.
+absorbleft = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mes
+nbregions = 3 # nb of regions and model number for each
+1 60 1 28 1
+30 30 26 28 2
+31 31 26 28 3
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/README
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/README (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/README 2011-12-29 00:40:50 UTC (rev 19324)
@@ -0,0 +1,29 @@
+----------------------------------------------------------------------
+README
+----------------------------------------------------------------------
+
+This example propagates an initial (incident) Rayleigh wave along the surface of a 2D homogeneous medium with a crack, similar to Figure 3 of Komatitsch and Tromp (1999) but with a crack instead of a canyon.
+
+TO RUN:
+
+0. Read the user manual in SPECFEM2D/doc/manual_SPECFEM2D.pdf
+
+1. in the SPECFEM2D root directory, configure, e.g.,
+ ./configure FC=gfortran
+
+2. compile:
+ make all
+
+3. cd EXAMPLES/Rayleigh_wave_with_crack
+
+4. execute script to run mesher and solver for the PSV case:
+ ./process.sh
+
+5. check out the output files in the local directory OUTPUT_FILES
+
+References:
+-----------
+
+Dimitri Komatitsch, Jean-Pierre Vilotte, Rossana Vai, José M. Castillo-Covarrubias and Francisco J. Sánchez-Sesma, The Spectral Element method for elastic wave equations: application to 2D and 3D seismic problems, International Journal for Numerical Methods in Engineering, vol. 45, p. 1139-1164 (1999).
+
+----------------------------------------------------------------------
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/SOURCE_Rayleigh_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/SOURCE_Rayleigh_2D (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/SOURCE_Rayleigh_2D 2011-12-29 00:40:50 UTC (rev 19324)
@@ -0,0 +1,13 @@
+# source parameters
+source_surf = .false. # source inside the medium or at the surface
+xs = 3.5 # source location x in meters
+zs = 9. # source location z in meters
+source_type = 3 # 1 for plane P waves, 2 for plane SV waves, 3 for Rayleigh wave
+time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0 = 1.0 # dominant source frequency (Hz) if not Dirac or Heaviside
+tshift = 0.0 # time shift when multi sources (if one source, must be zero)
+angleforce = 0. # angle of the source (for a force only)
+Mxx = 1. # Mxx component (for a moment tensor source only)
+Mzz = 1. # Mzz component (for a moment tensor source only)
+Mxz = 0. # Mxz component (for a moment tensor source only)
+factor = 1.d10 # amplification factor
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/interfaces_Rayleigh_flat.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/interfaces_Rayleigh_flat.dat (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/interfaces_Rayleigh_flat.dat 2011-12-29 00:40:50 UTC (rev 19324)
@@ -0,0 +1,26 @@
+#
+# number of interfaces
+#
+ 2
+#
+# for each interface below, we give the number of points and then x,y for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 19 0
+#
+# interface number 2 (topography, top of the mesh)
+#
+ 2
+ 0 9
+ 19 9
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+ 28
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/process.sh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/process.sh (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/process.sh 2011-12-29 00:40:50 UTC (rev 19324)
@@ -0,0 +1,69 @@
+#!/bin/bash
+#
+# script runs mesher and solver (in serial)
+# using this example setup
+#
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo
+echo "(will take a few minutes)"
+echo
+
+# sets up directory structure in current example directoy
+echo
+echo " setting up example..."
+echo
+
+mkdir -p OUTPUT_FILES
+mkdir -p DATA
+
+# sets up local DATA/ directory
+cd DATA/
+cp ../Par_file_Rayleigh_2D Par_file
+cp ../interfaces_Rayleigh_flat.dat .
+cp ../SOURCE_Rayleigh_2D SOURCE
+cd ../
+
+# cleans output files
+rm -rf OUTPUT_FILES/*
+
+# compiles executables in root directory
+cd ../../
+make > tmp.log
+cd $currentdir
+
+# links executables
+rm -f xmeshfem2D xspecfem2D
+ln -s ../../bin/xmeshfem2D
+ln -s ../../bin/xspecfem2D
+
+# stores setup
+cp DATA/Par_file OUTPUT_FILES/
+cp DATA/SOURCE OUTPUT_FILES/
+
+# runs database generation
+echo
+echo " running mesher..."
+echo
+./xmeshfem2D > OUTPUT_FILES/output_mesher.txt
+./xmeshfem2D
+
+# runs simulation
+echo
+echo " running solver..."
+echo
+#./xspecfem2D > OUTPUT_FILES/output_solver.txt
+./xspecfem2D
+
+# stores output
+cp DATA/SOURCE_xz.dat OUTPUT_FILES/
+cp DATA/STATIONS OUTPUT_FILES/
+cp DATA/STATIONS_target OUTPUT_FILES/
+
+echo
+echo "see results in directory: OUTPUT_FILES/"
+echo
+echo "done"
+echo `date`
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave_with_crack/process.sh
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_regions.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_regions.f90 2011-12-28 14:36:37 UTC (rev 19323)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_regions.f90 2011-12-29 00:40:50 UTC (rev 19324)
@@ -135,10 +135,10 @@
print *,' -----'
! store density and velocity model
- do i = ixdebregion,ixfinregion
- do j = izdebregion,izfinregion
- num_material((j-1)*nxread+i) = imaterial_number
- enddo
+ do j = izdebregion,izfinregion
+ do i = ixdebregion,ixfinregion
+ num_material((j-1)*nxread+i) = imaterial_number
+ enddo
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-12-28 14:36:37 UTC (rev 19323)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-12-29 00:40:50 UTC (rev 19324)
@@ -401,6 +401,13 @@
logical :: anyabs
double precision :: dxd,dyd,dzd,dcurld,valux,valuy,valuz,valcurl,hlagrange,rhol,xi,gamma,x,z
+!! DK DK Dec 2011: add a crack manually
+ logical, parameter :: ADD_A_CRACK = .true. ! .false.
+!! must be set equal to the number of spectral elements on one vertical side of the crack
+ integer :: NB_POINTS_TO_ADD_TO_NPGEO = 3
+ integer :: check_nb_points_to_add_to_npgeo,current_last_point,npgeo_ori,original_value,ignod
+ logical :: already_found_a_crack_element
+
! coefficients of the explicit Newmark time scheme
integer NSTEP
double precision :: deltatover2,deltatsquareover2,time
@@ -908,6 +915,9 @@
ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
if(nproc_read_from_database < 1) stop 'should have nproc_read_from_database >= 1'
if(nproc /= nproc_read_from_database) stop 'must always have nproc == nproc_read_from_database'
+!! DK DK Dec 2011: add a crack manually
+ npgeo_ori = npgeo
+ if(ADD_A_CRACK) npgeo = npgeo + NB_POINTS_TO_ADD_TO_NPGEO
#ifndef USE_MPI
if(PERFORM_CUTHILL_MCKEE) then
@@ -921,7 +931,8 @@
do ipass = 1,NUMBER_OF_PASSES
! starts reading in Database file
- if(ipass > 1) call read_databases_init(myrank,ipass, &
+ if(ipass > 1) &
+ call read_databases_init(myrank,ipass, &
simulation_title,SIMULATION_TYPE,NOISE_TOMOGRAPHY,SAVE_FORWARD,npgeo,nproc_read_from_database, &
gnuplot,interpol,NTSTEP_BETWEEN_OUTPUT_INFO, &
output_postscript_snapshot,output_color_image,colors,numbers, &
@@ -988,7 +999,9 @@
! reads the spectral macrobloc nodal coordinates
! and basic properties of the spectral elements
- call read_databases_coorg_elem(myrank,ipass,npgeo,coorg,numat,ngnod,nspec, &
+!! DK DK call read_databases_coorg_elem(myrank,ipass,npgeo,coorg,numat,ngnod,nspec, &
+!! DK DK Dec 2011: added a crack manually
+ call read_databases_coorg_elem(myrank,ipass,npgeo_ori,coorg,numat,ngnod,nspec, &
pointsdisp,plot_lowerleft_corner_only, &
nelemabs,nelem_acoustic_surface, &
num_fluid_solid_edges,num_fluid_poro_edges, &
@@ -1052,7 +1065,90 @@
call read_databases_mato(ipass,nspec,ngnod,kmato,knods, &
perm,antecedent_list)
+!! DK DK Dec 2011: handle a crack added manually
+ if(ADD_A_CRACK) then
+#ifdef USE_MPI
+ stop 'currently only serial runs are handled when adding a crack manually'
+#endif
+!! DK DK material number 2 indicates the spectral elements that form the left vertical side of the crack
+ check_nb_points_to_add_to_npgeo = count(kmato == 2)
+ print *
+ print *,'adding a crack manually'
+ print *,'need to add ',nb_points_to_add_to_npgeo,' npgeo mesh points to do that'
+
+ if(check_nb_points_to_add_to_npgeo /= NB_POINTS_TO_ADD_TO_NPGEO) &
+ stop 'must have check_nb_points_to_add_to_npgeo == NB_POINTS_TO_ADD_TO_NPGEO when adding a crack manually'
+
+ if(ngnod /= 4) stop 'must currently have ngnod == 4 when adding a crack manually'
+
+ if(PERFORM_CUTHILL_MCKEE) stop 'must not have PERFORM_CUTHILL_MCKEE when adding a crack manually'
+
+ if(FAST_NUMBERING) stop 'must not have FAST_NUMBERING when adding a crack manually'
+
+!! DK DK modify arrays "knods" and "coorg" to introduce the crack manually by duplicating and splitting the nodes
+ already_found_a_crack_element = .false.
+ current_last_point = npgeo_ori
+
+ do ispec = 1,nspec-1
+!! DK DK my convention is to introduce a vertical crack between two elements with material numbers 2 and 3
+ if(kmato(ispec) == 2 .and. kmato(ispec+1) == 3) then
+
+ print *,'adding a crack between elements ',ispec,' and ',ispec+1
+
+!! DK DK duplicate and split the lower-right corner of this element,
+!! DK DK except if it is the first crack element found, because then it is the crack
+!! DK DK tip and thus it should be assembled rather than split.
+!! DK DK Lower-right corner of an element is local npgeo point #2
+ if(already_found_a_crack_element .and. knods(2,ispec) <= npgeo_ori) then
+ current_last_point = current_last_point + 1
+ original_value = knods(2,ispec)
+!! DK DK split this point number in all the elements in which it appears
+ do ispec2 = 1,nspec
+! do this only for elements that define the left vertical edge of the crack
+ if(kmato(ispec2) /= 2) cycle
+ do ignod = 1,ngnod
+ if(knods(ignod,ispec2) == original_value) then
+ knods(ignod,ispec2) = current_last_point
+ coorg(:,current_last_point) = coorg(:,original_value)
+ endif
+ enddo
+ enddo
+ endif
+
+!! DK DK duplicate and split the upper-right corner of this element
+ already_found_a_crack_element = .true.
+
+!! DK DK Upper-right corner of an element is local npgeo point #3
+ if(knods(3,ispec) <= npgeo_ori) then
+
+ current_last_point = current_last_point + 1
+ original_value = knods(3,ispec)
+!! DK DK split this point number in all the elements in which it appears
+ do ispec2 = 1,nspec
+! do this only for elements that define the left vertical edge of the crack
+ if(kmato(ispec2) /= 2) cycle
+ do ignod = 1,ngnod
+ if(knods(ignod,ispec2) == original_value) then
+ knods(ignod,ispec2) = current_last_point
+ coorg(:,current_last_point) = coorg(:,original_value)
+ endif
+ enddo
+ enddo
+ endif
+
+ endif ! of if(kmato(ispec) == 2 .and. kmato(ispec+1) == 3)
+
+ enddo
+
+ if(current_last_point /= npgeo) then
+ print *,'current_last_point = ',current_last_point
+ print *,'npgeo_new = ',npgeo
+ stop 'did not find the right total number of points, should have current_last_point == npgeo_new'
+ endif
+
+ endif ! of if(ADD_A_CRACK) then
+
!-------------------------------------------------------------------------------
!---- determine if each spectral element is elastic, poroelastic, or acoustic
!-------------------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list