[cig-commits] r14285 - seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Mar 10 07:45:41 PDT 2009
Author: dkomati1
Date: 2009-03-10 07:45:41 -0700 (Tue, 10 Mar 2009)
New Revision: 14285
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90
Log:
added support for Piero Basini's new mesh of Antarctica
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 08:07:03 UTC (rev 14284)
+++ 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)
@@ -1,798 +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.
-!
-!=====================================================================
-
-! 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_refined.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