[cig-commits] r14134 - in seismo/3D/SPECFEM3D_SESAME/trunk: . UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Feb 24 06:50:07 PST 2009
Author: dkomati1
Date: 2009-02-24 06:50:07 -0800 (Tue, 24 Feb 2009)
New Revision: 14134
Added:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90
Removed:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/check_mesh_quality_AVS_DX.f90
seismo/3D/SPECFEM3D_SESAME/trunk/write_AVS_DX_mesh_quality_data.f90
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
Log:
suppressed all the old "mesh_quality_AVS_DX" routines
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90 2009-02-24 02:04:08 UTC (rev 14133)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90 2009-02-24 14:50:07 UTC (rev 14134)
@@ -1,111 +0,0 @@
-
- program analyze_CUBIT_Abaqus_mesh
-
-! Dimitri Komatitsch, University of Pau, France, February 2009.
-
-! Analyze a CUBIT mesh stored in Abaqus ASCII format and display the minimum
-! and maximum length of an edge.
-
- implicit none
-
-! number of points and of hex or quad elements
-! number of points of a hex or quad element
-
-! case of HOMOGENE_3D_in_meters.inp
-! integer, parameter :: NPOIN = 21212, NELEM = 18576, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-
-! case of regolite_3D_rego3d_70m_in_meters.inp
-! integer, parameter :: NPOIN = 4050696, NELEM = 3410265, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-
-! case of HOMOGENE_2D_in_meters.inp
-! integer, parameter :: NPOIN = 3882, NELEM = 3744, NGNOD = 4
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-
-! case of Eros_complet_regolite_fractures_2D_in_meters.inp
- integer, parameter :: NPOIN = 57630, NELEM = 107904, NGNOD = 4
- logical, parameter :: IGNORE_OTHER_HEADERS = .true.
-
- real, dimension(NPOIN) :: x,y,z
-
- integer, dimension(NGNOD,NELEM) :: ibool
-
- integer :: i,j,ielem,iread
-
- real :: distmin,distmax,dist
-
-! read the mesh
-
- print *,'start reading the CUBIT file'
-
-! skip the header
- read(*,*)
- read(*,*)
- read(*,*)
-
-! read the points
- do i = 1,NPOIN
- read(*,*) iread,x(i),y(i),z(i)
- if(iread /= i) then
- print *,'error at i,iread = ',i,iread
- stop 'wrong input for a point'
- endif
- enddo
-
-! skip the header
- read(*,*)
-
-! read the elements
- do i = 1,NELEM
-
- if(NGNOD == 4) then
-
-!! DK DK ignore other headers for 2D mesh of Eros with fractures, which has multiple material sets
- if(IGNORE_OTHER_HEADERS .and. (i == 51099 .or. i == 105264)) read(*,*)
-
- read(*,*) iread,ibool(1,i),ibool(2,i),ibool(3,i),ibool(4,i)
- else if(NGNOD == 8) then
- read(*,*) iread,ibool(1,i),ibool(2,i),ibool(3,i),ibool(4,i),ibool(5,i),ibool(6,i),ibool(7,i),ibool(8,i)
- else
- stop 'NGNOD should be 4 or 8'
- endif
-
- if(iread /= i) then
- print *,'error at i,iread = ',i,iread
- stop 'wrong input for an element'
- endif
-
- enddo
-
- print *,'end reading the CUBIT file'
- print *,'start computing the minimum and maximum edge size'
-
-! compute minimum and maximum distance
- distmin = + 1.e30
- distmax = - 1.e30
-
-! loop on all the elements
- do ielem = 1,NELEM
- if(mod(ielem,10000) == 0) print *,'processed ',ielem,' elements out of ',NELEM
- do i = 1,NGNOD
- do j = 1,NGNOD
- if(i /= j) then
- dist = sqrt((x(ibool(i,ielem)) - x(ibool(j,ielem)))**2 + &
- (y(ibool(i,ielem)) - y(ibool(j,ielem)))**2 + &
- (z(ibool(i,ielem)) - z(ibool(j,ielem)))**2)
- if(dist < distmin) distmin = dist
- if(dist > distmax) distmax = dist
- endif
- enddo
- enddo
- enddo
-
- print *
- print *,'minimum length of an edge in the whole mesh (m) = ',distmin
- print *,'minimum length divided in 4 average segments (m) = ',distmin/4.
- print *,'maximum length of an edge in the whole mesh (m) = ',distmax
- print *
-
- end program analyze_CUBIT_Abaqus_mesh
-
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 (from rev 14132, seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 2009-02-24 14:50:07 UTC (rev 14134)
@@ -0,0 +1,111 @@
+
+ program analyze_CUBIT_Abaqus_mesh
+
+! Dimitri Komatitsch, University of Pau, France, February 2009.
+
+! Analyze a CUBIT mesh stored in Abaqus ASCII format and display the minimum
+! and maximum length of an edge.
+
+ implicit none
+
+! number of points and of hex or quad elements
+! number of points of a hex or quad element
+
+! case of HOMOGENE_3D_in_meters.inp
+! integer, parameter :: NPOIN = 21212, NELEM = 18576, NGNOD = 8
+! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+
+! case of regolite_3D_rego3d_70m_in_meters.inp
+! integer, parameter :: NPOIN = 4050696, NELEM = 3410265, NGNOD = 8
+! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+
+! case of HOMOGENE_2D_in_meters.inp
+! integer, parameter :: NPOIN = 3882, NELEM = 3744, NGNOD = 4
+! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+
+! case of Eros_complet_regolite_fractures_2D_in_meters.inp
+ integer, parameter :: NPOIN = 57630, NELEM = 107904, NGNOD = 4
+ logical, parameter :: IGNORE_OTHER_HEADERS = .true.
+
+ real, dimension(NPOIN) :: x,y,z
+
+ integer, dimension(NGNOD,NELEM) :: ibool
+
+ integer :: i,j,ielem,iread
+
+ real :: distmin,distmax,dist
+
+! read the mesh
+
+ print *,'start reading the CUBIT file'
+
+! skip the header
+ read(*,*)
+ read(*,*)
+ read(*,*)
+
+! read the points
+ do i = 1,NPOIN
+ read(*,*) iread,x(i),y(i),z(i)
+ if(iread /= i) then
+ print *,'error at i,iread = ',i,iread
+ stop 'wrong input for a point'
+ endif
+ enddo
+
+! skip the header
+ read(*,*)
+
+! read the elements
+ do i = 1,NELEM
+
+ if(NGNOD == 4) then
+
+!! DK DK ignore other headers for 2D mesh of Eros with fractures, which has multiple material sets
+ if(IGNORE_OTHER_HEADERS .and. (i == 51099 .or. i == 105264)) read(*,*)
+
+ read(*,*) iread,ibool(1,i),ibool(2,i),ibool(3,i),ibool(4,i)
+ else if(NGNOD == 8) then
+ read(*,*) iread,ibool(1,i),ibool(2,i),ibool(3,i),ibool(4,i),ibool(5,i),ibool(6,i),ibool(7,i),ibool(8,i)
+ else
+ stop 'NGNOD should be 4 or 8'
+ endif
+
+ if(iread /= i) then
+ print *,'error at i,iread = ',i,iread
+ stop 'wrong input for an element'
+ endif
+
+ enddo
+
+ print *,'end reading the CUBIT file'
+ print *,'start computing the minimum and maximum edge size'
+
+! compute minimum and maximum distance
+ distmin = + 1.e30
+ distmax = - 1.e30
+
+! loop on all the elements
+ do ielem = 1,NELEM
+ if(mod(ielem,10000) == 0) print *,'processed ',ielem,' elements out of ',NELEM
+ do i = 1,NGNOD
+ do j = 1,NGNOD
+ if(i /= j) then
+ dist = sqrt((x(ibool(i,ielem)) - x(ibool(j,ielem)))**2 + &
+ (y(ibool(i,ielem)) - y(ibool(j,ielem)))**2 + &
+ (z(ibool(i,ielem)) - z(ibool(j,ielem)))**2)
+ if(dist < distmin) distmin = dist
+ if(dist > distmax) distmax = dist
+ endif
+ enddo
+ enddo
+ enddo
+
+ print *
+ print *,'minimum length of an edge in the whole mesh (m) = ',distmin
+ print *,'minimum length divided in 4 average segments (m) = ',distmin/4.
+ print *,'maximum length of an edge in the whole mesh (m) = ',distmax
+ print *
+
+ end program analyze_CUBIT_Abaqus_mesh
+
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90
___________________________________________________________________
Name: svn:mergeinfo
+
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/check_mesh_quality_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/check_mesh_quality_AVS_DX.f90 2009-02-24 02:04:08 UTC (rev 14133)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/check_mesh_quality_AVS_DX.f90 2009-02-24 14:50:07 UTC (rev 14134)
@@ -1,599 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! 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 2 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.
-!
-!=====================================================================
-
-! combine mesh quality data files to check the mesh
-! displays statistics on mesh quality
-! and creates an AVS or DX file showing a given range of elements
-
-!! DK DK
-!! DK DK this routine could be improved:
-!! DK DK
-!! DK DK - add full OpenDX support
-!! DK DK - debug aspect ratio
-!! DK DK - add mean in addition to min and max of ratios
-!! DK DK
-
- program mesh_quality_AVS_DX
-
- implicit none
-
- include "constants.h"
-
- integer iproc,nspec,npoin
- integer ispec
- integer iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
- integer ipoin,numpoin,iglobpointoffset,ntotpoin,ntotspec
- integer numelem,iglobelemoffset
- integer idoubling,iformat
- integer ntotpoinAVS_DX,ntotspecAVS_DX
-
- double precision xval,yval,zval
-
-! processor identification
- character(len=150) prname
-
-! parameters read from parameter file
- integer NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT, &
- NER_MOHO_16,NER_BOTTOM_MOHO,NEX_XI,NEX_ETA, &
- NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
- integer NSOURCES
-
- double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
- double precision DT,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,HDUR_MOVIE
- double precision THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
-
- logical HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,ATTENUATION,USE_OLSEN_ATTENUATION, &
- OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL, &
- BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS,SAVE_FORWARD
- logical ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
- logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
- USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
- integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
-
- character(len=150) OUTPUT_FILES,LOCAL_PATH,MODEL
-
-! parameters deduced from parameters read from file
- integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
- integer NER
-
-! for all the regions
- integer NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
- NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
-
-! for quality of mesh
- logical, dimension(:), allocatable :: mask_ibool
- double precision equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio
- double precision equiangle_skewness_min,edge_aspect_ratio_min,diagonal_aspect_ratio_min
- double precision equiangle_skewness_max,edge_aspect_ratio_max,diagonal_aspect_ratio_max
- double precision skewness_AVS_DX_min,skewness_AVS_DX_max
-
-! for stability and number of points per wavelength
- double precision stability,points_per_wavelength
- double precision stability_min,points_per_wavelength_min
- double precision stability_max,points_per_wavelength_max
-
-! for histogram
- integer, parameter :: NCLASS = 20
- integer classes_skewness(0:NCLASS-1)
- integer iclass
- double precision current_percent,total_percent
-
- integer proc_p1,proc_p2
-
-! ************** PROGRAM STARTS HERE **************
-
- print *
- print *,'Recombining all mesh quality files for slices'
- print *
-
- print *,'1 = create files in AVS UCD format'
- print *,'2 = create files in OpenDX format'
- print *,'any other value = exit'
- print *
- print *,'enter value:'
-!!!!!! DK DK read(5,*) iformat
-!! DK DK impose AVS format for now, OpenDX format not implemented yet
- iformat = 1
- if(iformat<1 .or. iformat>2) stop 'exiting...'
-
-! read range of skewness used for elements
- print *,'enter minimum skewness for AVS or DX (between 0. and 1.):'
- read(5,*) skewness_AVS_DX_min
- if(skewness_AVS_DX_min < 0.d0) skewness_AVS_DX_min = 0.d0
- if(skewness_AVS_DX_min > 0.99999d0) skewness_AVS_DX_min = 0.99999d0
-
- print *,'enter maximum skewness for AVS or DX (between 0. and 1.):'
- read(5,*) skewness_AVS_DX_max
- if(skewness_AVS_DX_max < 0.d0) skewness_AVS_DX_max = 0.d0
- if(skewness_AVS_DX_max > 0.99999d0) skewness_AVS_DX_max = 0.99999d0
-
- if(skewness_AVS_DX_min > skewness_AVS_DX_max) stop 'incorrect skewness range'
-
- print *
- print *,'reading parameter file'
- print *
-
-! read the parameter file
- call read_parameter_file(LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX, &
- UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
- NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
- NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,DT, &
- ATTENUATION,USE_OLSEN_ATTENUATION,HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,LOCAL_PATH,NSOURCES, &
- THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL,ANISOTROPY, &
- BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS, &
- MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
- NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
- SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD)
-
- if(.not. SAVE_MESH_FILES) stop 'AVS or DX files were not saved by the mesher'
-
-! compute other parameters based upon values read
- call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
- NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
- NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
- NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
-
-! get the base pathname for output files
- call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
-
- print *
- print *,'There are ',NPROC,' slices numbered from 0 to ',NPROC-1
- print *
-
-! use all the slices to determine correct statistics
- proc_p1 = 0
- proc_p2 = NPROC - 1
-
-! set total number of points and elements to zero
- ntotpoin = 0
- ntotspec = 0
-
-! loop on the selected range of processors
- do iproc = proc_p1,proc_p2
-
- print *,'Reading slice ',iproc
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,iproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXpoints.txt',status='old',action='read')
- read(10,*) npoin
- print *,'There are ',npoin,' global AVS or DX points in the slice'
- ntotpoin = ntotpoin + npoin
- close(10)
-
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelements.txt',status='old',action='read')
- read(10,*) nspec
- print *,'There are ',nspec,' AVS or DX elements in the slice'
- ntotspec = ntotspec + nspec
- close(10)
-
- enddo
-
- print *
- print *,'There is a total of ',ntotspec,' elements in all the slices'
- print *,'There is a total of ',ntotpoin,' points in all the slices'
- print *
-
-! ************* compute min and max of skewness and ratios ******************
-
-! erase minimum and maximum of quality numbers
- equiangle_skewness_min = + HUGEVAL
- edge_aspect_ratio_min = + HUGEVAL
- diagonal_aspect_ratio_min = + HUGEVAL
- stability_min = + HUGEVAL
- points_per_wavelength_min = + HUGEVAL
-
- equiangle_skewness_max = - HUGEVAL
- edge_aspect_ratio_max = - HUGEVAL
- diagonal_aspect_ratio_max = - HUGEVAL
- stability_max = - HUGEVAL
- points_per_wavelength_max = - HUGEVAL
-
-! set global element and point offsets to zero
- iglobelemoffset = 0
-
-! loop on the selected range of processors
- do iproc=proc_p1,proc_p2
-
- print *,'Reading slice ',iproc
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,iproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXmeshquality.txt',status='old',action='read')
-
- read(10,*) nspec
- print *,'There are ',nspec,' AVS or DX elements in the slice'
-
-! read local elements in this slice and output global AVS or DX elements
- do ispec=1,nspec
- read(10,*) numelem,equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,points_per_wavelength
- if(numelem /= ispec) stop 'incorrect element number'
-
-! mulitply stability number by time step
- stability = stability * DT
-
-! compute minimum and maximum of quality numbers
- equiangle_skewness_min = dmin1(equiangle_skewness_min,equiangle_skewness)
- edge_aspect_ratio_min = dmin1(edge_aspect_ratio_min,edge_aspect_ratio)
- diagonal_aspect_ratio_min = dmin1(diagonal_aspect_ratio_min,diagonal_aspect_ratio)
- stability_min = dmin1(stability_min,stability)
- points_per_wavelength_min = dmin1(points_per_wavelength_min,points_per_wavelength)
-
- equiangle_skewness_max = dmax1(equiangle_skewness_max,equiangle_skewness)
- edge_aspect_ratio_max = dmax1(edge_aspect_ratio_max,edge_aspect_ratio)
- diagonal_aspect_ratio_max = dmax1(diagonal_aspect_ratio_max,diagonal_aspect_ratio)
- stability_max = dmax1(stability_max,stability)
- points_per_wavelength_max = dmax1(points_per_wavelength_max,points_per_wavelength)
-
- enddo
-
- iglobelemoffset = iglobelemoffset + nspec
-
- close(10)
-
- enddo
-
- print *
- print *,'------------'
- print *,'mesh quality parameter definitions'
- print *
- print *,'equiangle skewness: 0. perfect 1. bad'
- print *,'skewness max deviation angle: 0. perfect 90. bad'
- print *,'skewness min mesh angle: 90. perfect 0. bad'
- print *,'edge aspect ratio: 1. perfect above 1. gives stretching factor'
- print *,'diagonal aspect ratio: 1. perfect above 1. gives stretching factor'
- print *,'------------'
-
- print *
- print *,'equiangle skewness max = ',equiangle_skewness_max
- print *,'equiangle skewness min = ',equiangle_skewness_min
- print *
- print *,'skewness max deviation angle = ',90.*equiangle_skewness_max
- print *,'skewness min mesh angle = ',90.*(1. - equiangle_skewness_max)
- print *
- print *,'edge aspect ratio max = ',edge_aspect_ratio_max
- print *,'edge aspect ratio min = ',edge_aspect_ratio_min
- print *
- print *,'diagonal aspect ratio max = ',diagonal_aspect_ratio_max
- print *,'diagonal aspect ratio min = ',diagonal_aspect_ratio_min
- print *
- print *,'stability max = ',stability_max
- print *,'stability min = ',stability_min
- print *
- print *,'points per wavelength max = ',points_per_wavelength_max
- print *,'points per wavelength min = ',points_per_wavelength_min
- print *
-
-! create statistics about mesh quality
-
- print *
- print *,'creating histogram and statistics of mesh quality - reading mesh data files'
- print *
-
-! erase histogram of skewness
- classes_skewness(:) = 0
-
-! erase number of elements belonging to skewness range for AVS_DX
- ntotspecAVS_DX = 0
-
-! set global element and point offsets to zero
- iglobelemoffset = 0
-
-! loop on the selected range of processors
- do iproc=proc_p1,proc_p2
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,iproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXmeshquality.txt',status='old',action='read')
-
- read(10,*) nspec
-
-! read local elements in this slice and output global AVS or DX elements
- do ispec=1,nspec
- read(10,*) numelem,equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,points_per_wavelength
- if(numelem /= ispec) stop 'incorrect element number'
-
-! store skewness in histogram
- iclass = int(equiangle_skewness * dble(NCLASS))
- if(iclass < 0) iclass = 0
- if(iclass > NCLASS-1) iclass = NCLASS-1
- classes_skewness(iclass) = classes_skewness(iclass) + 1
-
-! check if element belongs to requested skewness range
- if(equiangle_skewness >= skewness_AVS_DX_min .and. &
- equiangle_skewness <= skewness_AVS_DX_max) ntotspecAVS_DX = ntotspecAVS_DX + 1
-
- enddo
-
- iglobelemoffset = iglobelemoffset + nspec
-
- close(10)
-
- enddo
-
-! create histogram of skewness and save in Gnuplot file
- print *
- print *,'histogram of skewness (0. good - 1. bad):'
- print *
- total_percent = 0.
- open(unit=14,file=trim(OUTPUT_FILES)//'/mesh_quality_histogram.txt',status='unknown')
- do iclass = 0,NCLASS-1
- current_percent = 100.*dble(classes_skewness(iclass))/dble(ntotspec)
- total_percent = total_percent + current_percent
- print *,real(iclass/dble(NCLASS)),' - ',real((iclass+1)/dble(NCLASS)),classes_skewness(iclass),' ',sngl(current_percent),' %'
- write(14,*) 0.5*(real(iclass/dble(NCLASS)) + real((iclass+1)/dble(NCLASS))),' ',sngl(current_percent)
- enddo
- close(14)
-
-! create script for Gnuplot histogram file
- open(unit=14,file=trim(OUTPUT_FILES)//'/plot_mesh_quality_histogram.gnu',status='unknown')
- write(14,*) 'set term x11'
- write(14,*) 'set xrange [0:1]'
- write(14,*) 'set xtics 0,0.1,1'
- write(14,*) 'set boxwidth ',1./real(NCLASS)
- write(14,*) 'set xlabel "Skewness range"'
- write(14,*) 'set ylabel "Percentage of elements (%)"'
- write(14,*) 'plot "mesh_quality_histogram.txt" with boxes'
- write(14,*) 'pause -1 "hit any key..."'
- close(14)
-
- print *
- print *,'total number of elements = ',ntotspec
- print *,'total percentage = ',total_percent,' %'
- print *
-
-! display warning if maximum skewness is too high
- if(equiangle_skewness_max >= 0.80d0) then
- print *
- print *,'*********************************************'
- print *,'*********************************************'
- print *,' WARNING, mesh is bad (max skewness >= 0.80)'
- print *,'*********************************************'
- print *,'*********************************************'
- print *
- endif
-
-
-! ************* create AVS or DX file with elements in a certain range of skewness
-
- print *
- print *,'creating AVS or DX file with subset of elements in skewness range'
- print *,'between ',skewness_AVS_DX_min,' and ',skewness_AVS_DX_max
- print *
-
-! ************* count number of points without multiples ******************
-
-! set global element and point offsets to zero
- iglobpointoffset = 0
- iglobelemoffset = 0
-
-! allocate flag to remove multiples
- allocate(mask_ibool(ntotpoin))
- mask_ibool(:) = .false.
-
-! loop on the selected range of processors
- do iproc=proc_p1,proc_p2
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,iproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelements.txt',status='old',action='read')
- open(unit=12,file=prname(1:len_trim(prname))//'AVS_DXpoints.txt',status='old',action='read')
- open(unit=14,file=prname(1:len_trim(prname))//'AVS_DXmeshquality.txt',status='old',action='read')
-
- read(10,*) nspec
- read(12,*) npoin
- read(14,*) nspec
-
-! read local elements in this slice and output global AVS or DX elements
- do ispec=1,nspec
- read(10,*) numelem,idoubling,iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
- if(numelem /= ispec) stop 'incorrect element number'
-
- read(14,*) numelem,equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,points_per_wavelength
- if(numelem /= ispec) stop 'incorrect element number'
-
-! check if element belongs to requested skewness range
-! and flag all the points to remove multiples
- iglob1 = iglob1 + iglobpointoffset
- iglob2 = iglob2 + iglobpointoffset
- iglob3 = iglob3 + iglobpointoffset
- iglob4 = iglob4 + iglobpointoffset
- iglob5 = iglob5 + iglobpointoffset
- iglob6 = iglob6 + iglobpointoffset
- iglob7 = iglob7 + iglobpointoffset
- iglob8 = iglob8 + iglobpointoffset
- if(equiangle_skewness >= skewness_AVS_DX_min .and. equiangle_skewness <= skewness_AVS_DX_max) then
- mask_ibool(iglob1) = .true.
- mask_ibool(iglob2) = .true.
- mask_ibool(iglob3) = .true.
- mask_ibool(iglob4) = .true.
- mask_ibool(iglob5) = .true.
- mask_ibool(iglob6) = .true.
- mask_ibool(iglob7) = .true.
- mask_ibool(iglob8) = .true.
- endif
-
- enddo
-
- iglobelemoffset = iglobelemoffset + nspec
- iglobpointoffset = iglobpointoffset + npoin
-
- close(10)
- close(12)
- close(14)
-
- enddo
-
- if(ntotspecAVS_DX == 0) stop 'no elements in skewness range, no file created'
-
-! count number of independent points
- ntotpoinAVS_DX = count(mask_ibool(:))
-
- open(unit=11,file='AVS_meshquality.inp',status='unknown')
-
-! write AVS or DX header with element data
- write(11,*) ntotpoinAVS_DX,' ',ntotspecAVS_DX,' 0 1 0'
-
-! ************* generate points ******************
-
-! set global point offset to zero
- iglobpointoffset = 0
-
-! loop on the selected range of processors
- do iproc=proc_p1,proc_p2
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,iproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXpoints.txt',status='old',action='read')
- read(10,*) npoin
-
-! read local points in this slice and output global AVS or DX points
- do ipoin=1,npoin
- read(10,*) numpoin,xval,yval,zval
- if(numpoin /= ipoin) stop 'incorrect point number'
-! write to AVS or DX global file with correct offset if point has been selected
- if(mask_ibool(numpoin + iglobpointoffset)) &
- write(11,*) numpoin + iglobpointoffset,' ',sngl(xval),' ',sngl(yval),' ',sngl(zval)
-
- enddo
-
- iglobpointoffset = iglobpointoffset + npoin
-
- close(10)
-
- enddo
-
-! ************* generate elements ******************
-
-! set global element and point offsets to zero
- iglobpointoffset = 0
- iglobelemoffset = 0
-
-! loop on the selected range of processors
- do iproc=proc_p1,proc_p2
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,iproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelements.txt',status='old',action='read')
- open(unit=12,file=prname(1:len_trim(prname))//'AVS_DXpoints.txt',status='old',action='read')
- open(unit=14,file=prname(1:len_trim(prname))//'AVS_DXmeshquality.txt',status='old',action='read')
-
- read(10,*) nspec
- read(12,*) npoin
- read(14,*) nspec
-
-! read local elements in this slice and output global AVS or DX elements
- do ispec=1,nspec
- read(10,*) numelem,idoubling,iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
- if(numelem /= ispec) stop 'incorrect element number'
-
- read(14,*) numelem,equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,points_per_wavelength
- if(numelem /= ispec) stop 'incorrect element number'
-
-! write to AVS or DX global file with correct offset for hexahedra (3-D)
-! check if element belongs to requested skewness range
- iglob1 = iglob1 + iglobpointoffset
- iglob2 = iglob2 + iglobpointoffset
- iglob3 = iglob3 + iglobpointoffset
- iglob4 = iglob4 + iglobpointoffset
- iglob5 = iglob5 + iglobpointoffset
- iglob6 = iglob6 + iglobpointoffset
- iglob7 = iglob7 + iglobpointoffset
- iglob8 = iglob8 + iglobpointoffset
- if(equiangle_skewness >= skewness_AVS_DX_min .and. equiangle_skewness <= skewness_AVS_DX_max) &
- write(11,"(i6,' 1 hex ',i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6)") &
- numelem + iglobelemoffset,iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
- enddo
-
- iglobelemoffset = iglobelemoffset + nspec
- iglobpointoffset = iglobpointoffset + npoin
-
- close(10)
- close(12)
- close(14)
-
- enddo
-
-! ************* generate element data values ******************
-
-! output AVS or DX header for data
- write(11,*) '1 1'
- write(11,*) 'Zcoord, meters'
-
-! set global element and point offsets to zero
- iglobelemoffset = 0
-
-! loop on the selected range of processors
- do iproc=proc_p1,proc_p2
-
-! create the name for the database of the current slide
- call create_serial_name_database(prname,iproc,LOCAL_PATH,NPROC,OUTPUT_FILES)
-
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXmeshquality.txt',status='old',action='read')
-
- read(10,*) nspec
-
-! read local elements in this slice and output global AVS or DX elements
- do ispec=1,nspec
- read(10,*) numelem,equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,points_per_wavelength
- if(numelem /= ispec) stop 'incorrect element number'
-
-! write skewness data to AVS or DX global file with correct offset
-! scale skewness to [0:255] for AVS or DX color palette
-! check if element belongs to requested skewness range
- if(equiangle_skewness >= skewness_AVS_DX_min .and. equiangle_skewness <= skewness_AVS_DX_max) &
- write(11,*) numelem + iglobelemoffset,' ',255.*sngl(equiangle_skewness)
-
- enddo
-
- iglobelemoffset = iglobelemoffset + nspec
-
- close(10)
-
- enddo
-
-! close AVS or DX file
- close(11)
-
- print *
- print *,'there are ',ntotspecAVS_DX,' elements in AVS or DX skewness range ',skewness_AVS_DX_min,skewness_AVS_DX_max
- print *
-
- end program mesh_quality_AVS_DX
-
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-02-24 02:04:08 UTC (rev 14133)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-02-24 14:50:07 UTC (rev 14134)
@@ -977,9 +977,6 @@
! create AVS or DX mesh data for the slice, edges and faces
if(SAVE_MESH_FILES) then
call write_AVS_DX_global_data(myrank,prname,nspec,ibool,idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
-! It might be not necessary to save the ansotropic perturbation mesh data (c11store ...), so I just keep as it was
- call write_AVS_DX_mesh_quality_data(prname,nspec,xstore,ystore,zstore, &
- kappastore,mustore,rhostore)
call write_AVS_DX_global_faces_data(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool, &
idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
call write_AVS_DX_surface_data(myrank,prname,nspec,iboun,ibool, &
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/write_AVS_DX_mesh_quality_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/write_AVS_DX_mesh_quality_data.f90 2009-02-24 02:04:08 UTC (rev 14133)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/write_AVS_DX_mesh_quality_data.f90 2009-02-24 14:50:07 UTC (rev 14134)
@@ -1,261 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! 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 2 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.
-!
-!=====================================================================
-
-! create mesh quality data for the slice, to be recombined in postprocessing
-
- subroutine write_AVS_DX_mesh_quality_data(prname,nspec, &
- xstore,ystore,zstore,kappastore,mustore,rhostore)
-
- implicit none
-
- include "constants.h"
-
- integer nspec
- double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: kappastore,mustore,rhostore
-
- integer ispec
-
- double precision, dimension(8) :: xelm,yelm,zelm,vp,vs
-
- integer iface,icorner,jcorner
-
- double precision vectorA_x,vectorA_y,vectorA_z
- double precision vectorB_x,vectorB_y,vectorB_z
- double precision norm_A,norm_B,angle_vectors
- double precision distmin,distmax,dist,dist1,dist2,dist3,dist4
- double precision equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio
- double precision vs_min,vp_max
-
-! for stability and number of points per wavelength
- double precision highest_frequency_source,stability,points_per_wavelength
-
-! maximum polynomial degree for which we can compute the stability condition
- integer, parameter :: NGLL_MAX_STABILITY = 15
- double precision percent_GLL(NGLL_MAX_STABILITY)
-
-! topology of faces of cube for skewness
- integer faces_topo(6,6)
-
-! processor identification
- character(len=150) prname
-
-! data and element files identical to AVS_DXpoints.txt and AVS_DXelements.txt
-! created in regular AVS or DX routine, therefore not created again here
-
-! define percentage of smallest distance between GLL points for NGLL points
-! percentages were computed by calling the GLL points routine for each degree
- percent_GLL(2) = 100.d0
- percent_GLL(3) = 50.d0
- percent_GLL(4) = 27.639320225002102d0
- percent_GLL(5) = 17.267316464601141d0
- percent_GLL(6) = 11.747233803526763d0
- percent_GLL(7) = 8.4888051860716516d0
- percent_GLL(8) = 6.4129925745196719d0
- percent_GLL(9) = 5.0121002294269914d0
- percent_GLL(10) = 4.0233045916770571d0
- percent_GLL(11) = 3.2999284795970416d0
- percent_GLL(12) = 2.7550363888558858d0
- percent_GLL(13) = 2.3345076678918053d0
- percent_GLL(14) = 2.0032477366369594d0
- percent_GLL(15) = 1.7377036748080721d0
-
-! convert to real percentage
- percent_GLL(:) = percent_GLL(:) / 100.d0
-
-! check that the degree is not above the threshold for list of percentages
- if(NGLLX > NGLL_MAX_STABILITY) stop 'degree too high to compute stability value'
-
-! define topology of faces of cube for skewness
-
-! face 1
- faces_topo(1,1) = 1
- faces_topo(1,2) = 2
- faces_topo(1,3) = 6
- faces_topo(1,4) = 5
-
-! face 2
- faces_topo(2,1) = 2
- faces_topo(2,2) = 3
- faces_topo(2,3) = 7
- faces_topo(2,4) = 6
-
-! face 3
- faces_topo(3,1) = 4
- faces_topo(3,2) = 3
- faces_topo(3,3) = 7
- faces_topo(3,4) = 8
-
-! face 4
- faces_topo(4,1) = 1
- faces_topo(4,2) = 5
- faces_topo(4,3) = 8
- faces_topo(4,4) = 4
-
-! face 5
- faces_topo(5,1) = 1
- faces_topo(5,2) = 2
- faces_topo(5,3) = 3
- faces_topo(5,4) = 4
-
-! face 6
- faces_topo(6,1) = 5
- faces_topo(6,2) = 6
- faces_topo(6,3) = 7
- faces_topo(6,4) = 8
-
-! define wraparound for angles for skewness calculation
- faces_topo(:,5) = faces_topo(:,1)
- faces_topo(:,6) = faces_topo(:,2)
-
-! writing mesh quality data for each element
- open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXmeshquality.txt',status='unknown')
-
-! number of elements in AVS or DX file
- write(10,*) nspec
-
-! output global AVS or DX elements
- do ispec=1,nspec
-
-! define the coordinates of the 8 corners of the element
- xelm(1) = xstore(1,1,1,ispec)
- yelm(1) = ystore(1,1,1,ispec)
- zelm(1) = zstore(1,1,1,ispec)
- vs(1) = sqrt(mustore(1,1,1,ispec) / rhostore(1,1,1,ispec))
- vp(1) = sqrt(kappastore(1,1,1,ispec)/rhostore(1,1,1,ispec) + 4.d0*vs(1)*vs(1)/3.d0)
-
- xelm(2) = xstore(NGLLX,1,1,ispec)
- yelm(2) = ystore(NGLLX,1,1,ispec)
- zelm(2) = zstore(NGLLX,1,1,ispec)
- vs(2) = sqrt(mustore(NGLLX,1,1,ispec) / rhostore(NGLLX,1,1,ispec))
- vp(2) = sqrt(kappastore(NGLLX,1,1,ispec)/rhostore(NGLLX,1,1,ispec) + 4.d0*vs(2)*vs(2)/3.d0)
-
- xelm(3) = xstore(NGLLX,NGLLY,1,ispec)
- yelm(3) = ystore(NGLLX,NGLLY,1,ispec)
- zelm(3) = zstore(NGLLX,NGLLY,1,ispec)
- vs(3) = sqrt(mustore(NGLLX,NGLLY,1,ispec) / rhostore(NGLLX,NGLLY,1,ispec))
- vp(3) = sqrt(kappastore(NGLLX,NGLLY,1,ispec)/rhostore(NGLLX,NGLLY,1,ispec) + 4.d0*vs(3)*vs(3)/3.d0)
-
- xelm(4) = xstore(1,NGLLY,1,ispec)
- yelm(4) = ystore(1,NGLLY,1,ispec)
- zelm(4) = zstore(1,NGLLY,1,ispec)
- vs(4) = sqrt(mustore(1,NGLLY,1,ispec) / rhostore(1,NGLLY,1,ispec))
- vp(4) = sqrt(kappastore(1,NGLLY,1,ispec)/rhostore(1,NGLLY,1,ispec) + 4.d0*vs(4)*vs(4)/3.d0)
-
- xelm(5) = xstore(1,1,NGLLZ,ispec)
- yelm(5) = ystore(1,1,NGLLZ,ispec)
- zelm(5) = zstore(1,1,NGLLZ,ispec)
- vs(5) = sqrt(mustore(1,1,NGLLZ,ispec) / rhostore(1,1,NGLLZ,ispec))
- vp(5) = sqrt(kappastore(1,1,NGLLZ,ispec)/rhostore(1,1,NGLLZ,ispec) + 4.d0*vs(5)*vs(5)/3.d0)
-
- xelm(6) = xstore(NGLLX,1,NGLLZ,ispec)
- yelm(6) = ystore(NGLLX,1,NGLLZ,ispec)
- zelm(6) = zstore(NGLLX,1,NGLLZ,ispec)
- vs(6) = sqrt(mustore(NGLLX,1,NGLLZ,ispec) / rhostore(NGLLX,1,NGLLZ,ispec))
- vp(6) = sqrt(kappastore(NGLLX,1,NGLLZ,ispec)/rhostore(NGLLX,1,NGLLZ,ispec) + 4.d0*vs(6)*vs(6)/3.d0)
-
- xelm(7) = xstore(NGLLX,NGLLY,NGLLZ,ispec)
- yelm(7) = ystore(NGLLX,NGLLY,NGLLZ,ispec)
- zelm(7) = zstore(NGLLX,NGLLY,NGLLZ,ispec)
- vs(7) = sqrt(mustore(NGLLX,NGLLY,NGLLZ,ispec) / rhostore(NGLLX,NGLLY,NGLLZ,ispec))
- vp(7) = sqrt(kappastore(NGLLX,NGLLY,NGLLZ,ispec)/rhostore(NGLLX,NGLLY,NGLLZ,ispec) + 4.d0*vs(7)*vs(7)/3.d0)
-
- xelm(8) = xstore(1,NGLLY,NGLLZ,ispec)
- yelm(8) = ystore(1,NGLLY,NGLLZ,ispec)
- zelm(8) = zstore(1,NGLLY,NGLLZ,ispec)
- vs(8) = sqrt(mustore(1,NGLLY,NGLLZ,ispec) / rhostore(1,NGLLY,NGLLZ,ispec))
- vp(8) = sqrt(kappastore(1,NGLLY,NGLLZ,ispec)/rhostore(1,NGLLY,NGLLZ,ispec) + 4.d0*vs(8)*vs(8)/3.d0)
-
-! compute minimum Vs and maximum Vp
- vs_min = minval(vs(:))
- vp_max = maxval(vp(:))
-
-! compute equiangle skewness (as defined in Fluent/Gambit manual)
- equiangle_skewness = - HUGEVAL
- do iface = 1,6
- do icorner = 1,4
-
-! first vector of angle
- vectorA_x = xelm(faces_topo(iface,icorner)) - xelm(faces_topo(iface,icorner+1))
- vectorA_y = yelm(faces_topo(iface,icorner)) - yelm(faces_topo(iface,icorner+1))
- vectorA_z = zelm(faces_topo(iface,icorner)) - zelm(faces_topo(iface,icorner+1))
-
-! second vector of angle
- vectorB_x = xelm(faces_topo(iface,icorner+2)) - xelm(faces_topo(iface,icorner+1))
- vectorB_y = yelm(faces_topo(iface,icorner+2)) - yelm(faces_topo(iface,icorner+1))
- vectorB_z = zelm(faces_topo(iface,icorner+2)) - zelm(faces_topo(iface,icorner+1))
-
-! norm of vectors A and B
- norm_A = dsqrt(vectorA_x**2 + vectorA_y**2 + vectorA_z**2)
- norm_B = dsqrt(vectorB_x**2 + vectorB_y**2 + vectorB_z**2)
-
-! angle formed by the two vectors
- angle_vectors = dacos((vectorA_x*vectorB_x + vectorA_y*vectorB_y + vectorA_z*vectorB_z) / (norm_A * norm_B))
-
-! compute equiangle skewness
- equiangle_skewness = dmax1(equiangle_skewness,dabs(2.d0 * angle_vectors - PI) / PI)
-
- enddo
- enddo
-
-! compute edge aspect ratio using the 8 corners of the element
- distmin = + HUGEVAL
- distmax = - HUGEVAL
- do icorner = 1,8
- do jcorner = icorner + 1,8
- dist = dsqrt((xelm(jcorner) - xelm(icorner))**2 + (yelm(jcorner) - yelm(icorner))**2 + (zelm(jcorner) - zelm(icorner))**2)
- distmin = dmin1(distmin,dist)
- distmax = dmax1(distmax,dist)
- enddo
- enddo
- edge_aspect_ratio = distmax / distmin
-
-! compute stability value and number of points per shortest S wavelength
-! stability is multiplied by time step in check_mesh_quality_AVS_DX.f90
-!! DK DK incomplete formula for now, frequency of source hard wired
-!! DK DK should instead take from cmt_file, read and multiply
-!! DK DK directly later in check_mesh_quality.f90
- highest_frequency_source = 1.
- stability = vp_max / (distmin * percent_GLL(NGLLX))
- points_per_wavelength = NGLLX * vs_min / (distmax * highest_frequency_source)
-
-! compute diagonal aspect ratio
- dist1 = dsqrt((xelm(1) - xelm(7))**2 + (yelm(1) - yelm(7))**2 + (zelm(1) - zelm(7))**2)
- dist2 = dsqrt((xelm(2) - xelm(8))**2 + (yelm(2) - yelm(8))**2 + (zelm(2) - zelm(8))**2)
- dist3 = dsqrt((xelm(3) - xelm(5))**2 + (yelm(3) - yelm(5))**2 + (zelm(3) - zelm(5))**2)
- dist4 = dsqrt((xelm(4) - xelm(6))**2 + (yelm(4) - yelm(6))**2 + (zelm(4) - zelm(6))**2)
- distmin = dmin1(dist1,dist2,dist3,dist4)
- distmax = dmax1(dist1,dist2,dist3,dist4)
- diagonal_aspect_ratio = distmax / distmin
-
-! write mesh quality information for each element
- write(10,*) ispec,equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,points_per_wavelength
-
- enddo
-
- close(10)
-
- end subroutine write_AVS_DX_mesh_quality_data
-
More information about the CIG-COMMITS
mailing list