[cig-commits] r14286 - seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Mar 10 08:00:00 PDT 2009
Author: dkomati1
Date: 2009-03-10 08:00:00 -0700 (Tue, 10 Mar 2009)
New Revision: 14286
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90
Log:
re-added the file because I realize that I erased it a few minutes ago
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 2009-03-10 14:45:41 UTC (rev 14285)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 2009-03-10 15:00:00 UTC (rev 14286)
@@ -0,0 +1,798 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+! read a 2D or 3D CUBIT mesh file and display statistics about mesh quality;
+! and create an OpenDX file showing a given range of elements or a single element
+
+! Dimitri Komatitsch, University of Pau, France, March 2009.
+
+!! DK DK
+!! DK DK this routine could be improved by computing the mean in addition to min and max of ratios
+!! DK DK
+
+ program check_mesh_quality_CUBIT_Abaqus
+
+ implicit none
+
+ include "constants.h"
+
+! number of points and of hex or quad elements
+! number of points of a hex or quad element
+
+! character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_3D_lisse_300_in_meters.inp'
+! integer, parameter :: NPOIN = 98692, NSPEC = 90585, NGNOD = 8
+! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+! double precision, parameter :: delta_t = 1.d-3
+! double precision, parameter :: VP_MAX = 3000.d0
+
+! character(len=100), parameter :: cubit_mesh_file = 'regolite_3D_rego3d_70m_in_meters.inp'
+! integer, parameter :: NPOIN = 4050696, NSPEC = 3410265, NGNOD = 8
+! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+! double precision, parameter :: delta_t = 3.d-4
+! double precision, parameter :: VP_MAX = 3000.d0
+
+! character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_2D_in_meters.inp'
+! integer, parameter :: NPOIN = 3882, NSPEC = 3744, NGNOD = 4
+! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+! double precision, parameter :: delta_t = 5.d-3
+! double precision, parameter :: VP_MAX = 3000.d0
+
+! character(len=100), parameter :: cubit_mesh_file = 'eros_complexe_2d_regolite_fractures_modifie_in_meters.inp'
+! integer, parameter :: NPOIN = 57807, NSPEC = 56983, NGNOD = 4
+! logical, parameter :: IGNORE_OTHER_HEADERS = .true.
+! double precision, parameter :: delta_t = 1.5d-4
+! double precision, parameter :: VP_MAX = 3000.d0
+
+! character(len=100), parameter :: cubit_mesh_file = 'REGOLITE_only_no_fractures_2D_in_meters.inp'
+! integer, parameter :: NPOIN = 32536, NSPEC = 31695, NGNOD = 4
+! logical, parameter :: IGNORE_OTHER_HEADERS = .true.
+! double precision, parameter :: delta_t = 1.5d-4
+! double precision, parameter :: VP_MAX = 900.d0 ! because the smallest element is in the regolith layer, not in the bedrock
+
+ character(len=100), parameter :: cubit_mesh_file = 'david_mesh_doubl_500_900_6layers.inp'
+ integer, parameter :: NPOIN = 4513255, NSPEC = 4379190, NGNOD = 8
+ logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+ double precision, parameter :: delta_t = 1.5d-4
+ double precision, parameter :: VP_MAX = 3000.d0
+
+ double precision, dimension(NPOIN) :: x,y,z
+
+ integer, dimension(NGNOD,NSPEC) :: ibool
+
+ integer :: i,ispec,iread,iformat,ispec_min_edge_length,ispec_max_edge_length, &
+ ispec_begin,ispec_end,ispec_to_output
+
+! for quality of mesh
+ 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,distance_min,distance_max
+ double precision :: distmin,distmax
+
+! for stability
+ double precision :: stability,stability_min,stability_max
+
+! for histogram
+ integer, parameter :: NCLASS = 20
+ integer classes_skewness(0:NCLASS-1)
+ integer :: iclass
+ double precision :: current_percent,total_percent
+
+! to export elements that have a certain skewness range to OpenDX
+ integer :: ntotspecAVS_DX
+ logical :: USE_OPENDX
+
+ if(NGNOD /= 4 .and. NGNOD /= 8) stop 'NGNOD should be 4 or 8'
+
+ print *
+ print *,'1 = output elements above a certain skewness threshold in OpenDX format'
+ print *,'2 = output a given element in OpenDX format'
+ print *,'3 = do not output any OpenDX file'
+ print *
+ print *,'enter value:'
+ read(5,*) iformat
+
+ if(iformat < 1 .or. iformat > 3) stop 'exiting...'
+
+ if(iformat == 1 .or. iformat == 2) then
+ USE_OPENDX = .true.
+ else
+ USE_OPENDX = .false.
+ endif
+
+ if(USE_OPENDX) then
+
+ if(iformat == 1) then
+
+! read range of skewness used for elements
+ print *,'enter minimum skewness for OpenDX (between 0. and 0.99):'
+ 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 OpenDX (between 0. and 1.):'
+!!!!!!!!!!!!! read(5,*) skewness_AVS_DX_max
+ skewness_AVS_DX_max = 0.99999d0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ 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'
+
+ else
+ print *,'enter the element number to output in OpenDX format between 1 and ',NSPEC
+ read(5,*) ispec_to_output
+ if(ispec_to_output < 1 .or. ispec_to_output > NSPEC) stop 'incorrect element number to output'
+ endif
+
+ endif
+
+! read the mesh
+ print *
+ print *,'start reading the CUBIT file: ',cubit_mesh_file(1:len_trim(cubit_mesh_file))
+ open(unit=10,file=cubit_mesh_file,status='unknown',action='read')
+
+! skip the header
+ read(10,*)
+ read(10,*)
+ read(10,*)
+
+! read the points
+ do i = 1,NPOIN
+ read(10,*) 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(10,*)
+
+! read the elements
+ do i = 1,NSPEC
+
+ 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. cubit_mesh_file == 'eros_complexe_2d_regolite_fractures_modifie_in_meters.inp' &
+ .and. i == 5709) read(10,*)
+
+ if(IGNORE_OTHER_HEADERS .and. cubit_mesh_file == 'REGOLITE_only_no_fractures_2D_in_meters.inp' &
+ .and. i == 28429) read(10,*)
+
+ read(10,*) iread,ibool(1,i),ibool(2,i),ibool(3,i),ibool(4,i)
+ else if(NGNOD == 8) then
+ read(10,*) 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)
+ endif
+
+ if(iread /= i) then
+ print *,'error at i,iread = ',i,iread
+ stop 'wrong input for an element'
+ endif
+
+ enddo
+ close(10)
+
+ print *,'done reading the CUBIT file'
+ print *
+
+ print *,'start computing the minimum and maximum edge size'
+
+! ************* 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
+ distance_min = + HUGEVAL
+
+ equiangle_skewness_max = - HUGEVAL
+ edge_aspect_ratio_max = - HUGEVAL
+ diagonal_aspect_ratio_max = - HUGEVAL
+ stability_max = - HUGEVAL
+ distance_max = - HUGEVAL
+
+ ispec_min_edge_length = -1
+ ispec_max_edge_length = -1
+
+! loop on all the elements
+ do ispec = 1,NSPEC
+
+ if(mod(ispec,100000) == 0) print *,'processed ',ispec,' elements out of ',NSPEC
+
+ if(NGNOD == 4) then
+ call create_mesh_quality_data_2D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ else
+ call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ endif
+
+! store element number in which the edge of minimum or maximum length is located
+ if(distmin < distance_min) ispec_min_edge_length = ispec
+ if(distmax > distance_max) ispec_max_edge_length = ispec
+
+! compute minimum and maximum of quality numbers
+ equiangle_skewness_min = min(equiangle_skewness_min,equiangle_skewness)
+ edge_aspect_ratio_min = min(edge_aspect_ratio_min,edge_aspect_ratio)
+ diagonal_aspect_ratio_min = min(diagonal_aspect_ratio_min,diagonal_aspect_ratio)
+ stability_min = min(stability_min,stability)
+ distance_min = min(distance_min,distmin)
+
+ equiangle_skewness_max = max(equiangle_skewness_max,equiangle_skewness)
+ edge_aspect_ratio_max = max(edge_aspect_ratio_max,edge_aspect_ratio)
+ diagonal_aspect_ratio_max = max(diagonal_aspect_ratio_max,diagonal_aspect_ratio)
+ stability_max = max(stability_max,stability)
+ distance_max = max(distance_max,distmax)
+
+ enddo
+ print *,'done processing ',NSPEC,' elements out of ',NSPEC
+
+ 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 *,'edge aspect ratio: 1. perfect above 1. gives stretching factor'
+ print *,'diagonal aspect ratio: 1. perfect above 1. gives stretching factor'
+ print *,'------------'
+
+ print *
+ print *,'minimum length of an edge in the whole mesh (m) = ',distance_min,' in element ',ispec_min_edge_length
+ print *
+ print *,'maximum length of an edge in the whole mesh (m) = ',distance_max,' in element ',ispec_max_edge_length
+ print *
+ print *,'max equiangle skewness = ',equiangle_skewness_max
+! print *,'min equiangle skewness = ',equiangle_skewness_min
+ print *
+ print *,'max deviation angle from a right angle (90 degrees) is therefore = ',90.*equiangle_skewness_max
+ print *
+ print *,'worst angle in the mesh is therefore ',90.*(1. - equiangle_skewness_max)
+ print *,'or ',180. - 90.*(1. - equiangle_skewness_max),' degrees'
+ print *
+ print *,'max edge aspect ratio = ',edge_aspect_ratio_max
+! print *,'min edge aspect ratio = ',edge_aspect_ratio_min
+ print *
+ print *,'max diagonal aspect ratio = ',diagonal_aspect_ratio_max
+! print *,'min diagonal aspect ratio = ',diagonal_aspect_ratio_min
+ print *
+ print *,'max stability = ',stability_max
+! print *,'min stability = ',stability_min
+
+! create statistics about mesh quality
+
+ print *
+ print *,'creating histogram and statistics of mesh quality'
+
+! erase histogram of skewness
+ classes_skewness(:) = 0
+
+! loop on all the elements
+ do ispec = 1,NSPEC
+
+ if(NGNOD == 4) then
+ call create_mesh_quality_data_2D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ else
+ call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ endif
+
+! 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
+
+ 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='mesh_quality_histogram.txt',status='unknown')
+ do iclass = 0,NCLASS-1
+ current_percent = 100.*dble(classes_skewness(iclass))/dble(NSPEC)
+ 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='plot_mesh_quality_histogram.gnu',status='unknown')
+ write(14,*) 'set term x11'
+ write(14,*) '#set term gif'
+ write(14,*) '#set output "mesh_quality_histogram.gif"'
+ write(14,*)
+ 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 = ',NSPEC
+ print *
+
+! display warning if maximum skewness is too high
+ if(equiangle_skewness_max >= 0.75d0) then
+ print *
+ print *,'*********************************************'
+ print *,'*********************************************'
+ print *,' WARNING, mesh is bad (max skewness >= 0.75)'
+ print *,'*********************************************'
+ print *,'*********************************************'
+ print *
+ endif
+
+ if(total_percent < 99.9d0 .or. total_percent > 100.1d0) then
+ print *,'total percentage = ',total_percent,' %'
+ stop 'total percentage should be 100%'
+ endif
+
+! ************* create OpenDX file with elements in a certain range of skewness
+
+ if(USE_OPENDX) then
+
+ print *
+ if(iformat == 1) then
+ print *,'creating OpenDX file with subset of elements in skewness range'
+ print *,'between ',skewness_AVS_DX_min,' and ',skewness_AVS_DX_max
+ else
+ print *,'creating OpenDX file with element #',ispec_to_output
+ endif
+ print *
+
+! ************* count number of elements in skewness range *************
+
+! erase number of elements belonging to skewness range for AVS_DX
+ ntotspecAVS_DX = 0
+
+! loop on all the elements
+ if(iformat == 1) then
+
+ do ispec = 1,NSPEC
+
+ if(NGNOD == 4) then
+ call create_mesh_quality_data_2D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ else
+ call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ endif
+
+! 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
+
+ else
+! outputing a single element
+ ntotspecAVS_DX = 1
+ endif
+
+ if(ntotspecAVS_DX == 0) then
+ stop 'no elements in skewness range, no file created'
+ else if(iformat == 1) then
+ print *
+ print *,'there are ',ntotspecAVS_DX,' elements in AVS or DX skewness range ',skewness_AVS_DX_min,skewness_AVS_DX_max
+ print *
+ endif
+
+ open(unit=11,file='DX_mesh_quality.dx',status='unknown')
+
+! ************* generate points ******************
+
+! write OpenDX header
+ write(11,*) 'object 1 class array type float rank 1 shape 3 items ',NPOIN,' data follows'
+
+! write all the points
+ do i = 1,NPOIN
+ write(11,*) sngl(x(i)),sngl(y(i)),sngl(z(i))
+ enddo
+
+! ************* generate elements ******************
+
+ write(11,*) 'object 2 class array type int rank 1 shape ',NGNOD,' items ',ntotspecAVS_DX,' data follows'
+
+! loop on all the elements
+ if(iformat == 1) then
+ ispec_begin = 1
+ ispec_end = NSPEC
+ else
+ ispec_begin = ispec_to_output
+ ispec_end = ispec_to_output
+ endif
+
+ do ispec = ispec_begin,ispec_end
+
+ if(NGNOD == 4) then
+ call create_mesh_quality_data_2D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ else
+ call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ endif
+
+! check if element needs to be output
+ if(iformat == 2 .or. (iformat == 1 .and. &
+ equiangle_skewness >= skewness_AVS_DX_min .and. equiangle_skewness <= skewness_AVS_DX_max)) then
+! 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
+ if(NGNOD == 4) then
+ write(11,"(i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9,1x,i9)") &
+ ibool(1,ispec)-1, ibool(4,ispec)-1, ibool(2,ispec)-1, ibool(3,ispec)-1
+ else
+ 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
+ endif
+ if(iformat == 1) print *,'element ',ispec,' belongs to the range and has skewness = ',sngl(equiangle_skewness)
+ endif
+
+ enddo
+
+! ************* generate element data values ******************
+
+! output OpenDX header for data
+ if(NGNOD == 4) then
+ write(11,*) 'attribute "element type" string "quads"'
+ else
+ write(11,*) 'attribute "element type" string "cubes"'
+ endif
+ write(11,*) 'attribute "ref" string "positions"'
+ write(11,*) 'object 3 class array type float rank 0 items ',ntotspecAVS_DX,' data follows'
+
+! loop on all the elements
+ do ispec = ispec_begin,ispec_end
+
+ if(NGNOD == 4) then
+ call create_mesh_quality_data_2D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ else
+ call create_mesh_quality_data_3D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+ endif
+
+! check if element needs to be output
+ if(iformat == 2 .or. (iformat == 1 .and. &
+ equiangle_skewness >= skewness_AVS_DX_min .and. equiangle_skewness <= skewness_AVS_DX_max)) &
+ write(11,*) sngl(equiangle_skewness)
+
+ 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 OpenDX file
+ close(11)
+
+ endif
+
+ end program check_mesh_quality_CUBIT_Abaqus
+
+!
+!=====================================================================
+!
+
+! create mesh quality data for a given 3D spectral element
+
+ subroutine create_mesh_quality_data_3D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: iface,icorner,ispec,NSPEC,NPOIN,NGNOD,i
+
+ double precision, dimension(NPOIN) :: x,y,z
+
+ integer, dimension(NGNOD,NSPEC) :: ibool
+
+ double precision, dimension(NGNOD) :: xelm,yelm,zelm
+
+ 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
+
+! for stability
+ double precision :: stability,VP_MAX,delta_t
+
+! maximum polynomial degree for which we can compute the stability condition
+ integer, parameter :: NGLL_MAX_STABILITY = 15
+ double precision, dimension(NGLL_MAX_STABILITY) :: percent_GLL
+
+! topology of faces of cube for skewness
+ integer faces_topo(6,6)
+
+! store the corners of this element for the skewness routine
+ do i = 1,NGNOD
+ xelm(i) = x(ibool(i,ispec))
+ yelm(i) = y(ibool(i,ispec))
+ zelm(i) = z(ibool(i,ispec))
+ enddo
+
+! 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)
+
+! compute equiangle skewness (as defined in Fluent/Gambit manual)
+! and compute edge aspect ratio using the corners of the element
+ distmin = + HUGEVAL
+ distmax = - HUGEVAL
+ 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 = sqrt(vectorA_x**2 + vectorA_y**2 + vectorA_z**2)
+ norm_B = sqrt(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 = max(equiangle_skewness,dabs(2.d0 * angle_vectors - PI) / PI)
+
+! compute min and max size of an edge
+ dist = sqrt(vectorA_x**2 + vectorA_y**2 + vectorA_z**2)
+
+ distmin = min(distmin,dist)
+ distmax = max(distmax,dist)
+
+ enddo
+ enddo
+
+! compute edge aspect ratio
+ edge_aspect_ratio = distmax / distmin
+
+ stability = delta_t * VP_MAX / (distmin * percent_GLL(NGLLX))
+
+! compute diagonal aspect ratio
+ dist1 = sqrt((xelm(1) - xelm(7))**2 + (yelm(1) - yelm(7))**2 + (zelm(1) - zelm(7))**2)
+ dist2 = sqrt((xelm(2) - xelm(8))**2 + (yelm(2) - yelm(8))**2 + (zelm(2) - zelm(8))**2)
+ dist3 = sqrt((xelm(3) - xelm(5))**2 + (yelm(3) - yelm(5))**2 + (zelm(3) - zelm(5))**2)
+ dist4 = sqrt((xelm(4) - xelm(6))**2 + (yelm(4) - yelm(6))**2 + (zelm(4) - zelm(6))**2)
+ diagonal_aspect_ratio = max(dist1,dist2,dist3,dist4) / min(dist1,dist2,dist3,dist4)
+
+ end subroutine create_mesh_quality_data_3D
+
+!
+!=====================================================================
+!
+
+! create mesh quality data for a given 2D spectral element
+
+ subroutine create_mesh_quality_data_2D(x,y,z,ibool,ispec,NSPEC,NPOIN,NGNOD,VP_MAX,delta_t, &
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: icorner,ispec,NSPEC,NPOIN,NGNOD,i
+
+ double precision, dimension(NPOIN) :: x,y,z
+
+ integer, dimension(NGNOD,NSPEC) :: ibool
+
+ double precision, dimension(NGNOD) :: xelm,yelm,zelm
+
+ 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
+ double precision equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio
+
+! for stability
+ double precision :: stability,VP_MAX,delta_t
+
+! maximum polynomial degree for which we can compute the stability condition
+ integer, parameter :: NGLL_MAX_STABILITY = 15
+ double precision, dimension(NGLL_MAX_STABILITY) :: percent_GLL
+
+! topology of faces of cube for skewness
+! only one face in 2D
+ integer faces_topo(6)
+
+! store the corners of this element for the skewness routine
+ do i = 1,NGNOD
+ xelm(i) = x(ibool(i,ispec))
+ yelm(i) = y(ibool(i,ispec))
+ zelm(i) = z(ibool(i,ispec))
+ enddo
+
+! 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
+
+! only one face in 2D
+ faces_topo(1) = 1
+ faces_topo(2) = 2
+ faces_topo(3) = 3
+ faces_topo(4) = 4
+
+! define wraparound for angles for skewness calculation
+ faces_topo(5) = faces_topo(1)
+ faces_topo(6) = faces_topo(2)
+
+! compute equiangle skewness (as defined in Fluent/Gambit manual)
+! and compute edge aspect ratio using the corners of the element
+ distmin = + HUGEVAL
+ distmax = - HUGEVAL
+ equiangle_skewness = - HUGEVAL
+
+ do icorner = 1,4
+
+! first vector of angle
+ vectorA_x = xelm(faces_topo(icorner)) - xelm(faces_topo(icorner+1))
+ vectorA_y = yelm(faces_topo(icorner)) - yelm(faces_topo(icorner+1))
+ vectorA_z = zelm(faces_topo(icorner)) - zelm(faces_topo(icorner+1))
+
+! second vector of angle
+ vectorB_x = xelm(faces_topo(icorner+2)) - xelm(faces_topo(icorner+1))
+ vectorB_y = yelm(faces_topo(icorner+2)) - yelm(faces_topo(icorner+1))
+ vectorB_z = zelm(faces_topo(icorner+2)) - zelm(faces_topo(icorner+1))
+
+! norm of vectors A and B
+ norm_A = sqrt(vectorA_x**2 + vectorA_y**2 + vectorA_z**2)
+ norm_B = sqrt(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 = max(equiangle_skewness,dabs(2.d0 * angle_vectors - PI) / PI)
+
+! compute min and max size of an edge
+ dist = sqrt(vectorA_x**2 + vectorA_y**2 + vectorA_z**2)
+
+ distmin = min(distmin,dist)
+ distmax = max(distmax,dist)
+
+ enddo
+
+! compute edge aspect ratio
+ edge_aspect_ratio = distmax / distmin
+
+ stability = delta_t * VP_MAX / (distmin * percent_GLL(NGLLX))
+
+! compute diagonal aspect ratio
+ dist1 = sqrt((xelm(1) - xelm(3))**2 + (yelm(1) - yelm(3))**2 + (zelm(1) - zelm(3))**2)
+ dist2 = sqrt((xelm(2) - xelm(4))**2 + (yelm(2) - yelm(4))**2 + (zelm(2) - zelm(4))**2)
+ diagonal_aspect_ratio = max(dist1,dist2) / min(dist1,dist2)
+
+ end subroutine create_mesh_quality_data_2D
+
More information about the CIG-COMMITS
mailing list