[cig-commits] r14135 - seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Feb 24 06:51:55 PST 2009
Author: dkomati1
Date: 2009-02-24 06:51:55 -0800 (Tue, 24 Feb 2009)
New Revision: 14135
Added:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/constants.h
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/multiply_CUBIT_Abaqus_mesh_by_1000.f90
Log:
added code to analyze in detail the quality of CUBIT/Abaqus meshes, based on my older "mesh_quality_AVS_DX" routines
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-02-24 14:50:07 UTC (rev 14134)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 2009-02-24 14:51:55 UTC (rev 14135)
@@ -1,52 +1,144 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
- program analyze_CUBIT_Abaqus_mesh
+! read a 2D or 3D CUBIT mesh file and display statistics about mesh quality;
+! and create an OpenDX file showing a given range of bad elements
! 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.
+!! DK DK
+!! DK DK this routine could be improved:
+!! DK DK
+!! DK DK - debug aspect ratio
+!! DK DK - add 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
-! case of HOMOGENE_3D_in_meters.inp
-! integer, parameter :: NPOIN = 21212, NELEM = 18576, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+ 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
-! case of regolite_3D_rego3d_70m_in_meters.inp
-! integer, parameter :: NPOIN = 4050696, NELEM = 3410265, NGNOD = 8
+! 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
-! case of HOMOGENE_2D_in_meters.inp
-! integer, parameter :: NPOIN = 3882, NELEM = 3744, NGNOD = 4
+! 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
-! case of Eros_complet_regolite_fractures_2D_in_meters.inp
- integer, parameter :: NPOIN = 57630, NELEM = 107904, NGNOD = 4
- logical, parameter :: IGNORE_OTHER_HEADERS = .true.
+! character(len=100), parameter :: cubit_mesh_file = 'Eros_complet_regolite_fractures_2D_in_meters.inp'
+! integer, parameter :: NPOIN = 57630, NSPEC = 107904, NGNOD = 4
+! logical, parameter :: IGNORE_OTHER_HEADERS = .true.
+! double precision, parameter :: delta_t = 1.5d-4
- real, dimension(NPOIN) :: x,y,z
+ double precision, parameter :: VP_MAX = 3000.d0
- integer, dimension(NGNOD,NELEM) :: ibool
+ double precision, dimension(NPOIN) :: x,y,z
- integer :: i,j,ielem,iread
+ integer, dimension(NGNOD,NSPEC) :: ibool
- real :: distmin,distmax,dist
+ integer :: i,j,ispec,iread,iformat
+ double precision :: distmin,distmax,dist
+
+! 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
+
+! for stability
+ double precision stability
+ double precision stability_min
+ double precision 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 *,'1 = output bad elements in OpenDX format'
+ print *,'2 = do not output bad elements'
+ print *
+ print *,'enter value:'
+ read(5,*) iformat
+ if(iformat<1 .or. iformat>2) stop 'exiting...'
+ if(iformat == 1) then
+ USE_OPENDX = .true.
+ else
+ USE_OPENDX = .false.
+ endif
+
+ if(USE_OPENDX) 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'
+
+ endif
+
! read the mesh
+ 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')
- print *,'start reading the CUBIT file'
-
! skip the header
- read(*,*)
- read(*,*)
- read(*,*)
+ read(10,*)
+ read(10,*)
+ read(10,*)
! read the points
do i = 1,NPOIN
- read(*,*) iread,x(i),y(i),z(i)
+ 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'
@@ -54,21 +146,19 @@
enddo
! skip the header
- read(*,*)
+ read(10,*)
! read the elements
- do i = 1,NELEM
+ 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. (i == 51099 .or. i == 105264)) read(*,*)
+ if(IGNORE_OTHER_HEADERS .and. (i == 51099 .or. i == 105264)) read(10,*)
- read(*,*) iread,ibool(1,i),ibool(2,i),ibool(3,i),ibool(4,i)
+ read(10,*) 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'
+ 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
@@ -77,6 +167,7 @@
endif
enddo
+ close(10)
print *,'end reading the CUBIT file'
print *,'start computing the minimum and maximum edge size'
@@ -86,26 +177,584 @@
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 ispec = 1,NSPEC
+ if(mod(ispec,10000) == 0) print *,'processed ',ispec,' elements out of ',NSPEC
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)
+ dist = sqrt((x(ibool(i,ispec)) - x(ibool(j,ispec)))**2 + &
+ (y(ibool(i,ispec)) - y(ibool(j,ispec)))**2 + &
+ (z(ibool(i,ispec)) - z(ibool(j,ispec)))**2)
if(dist < distmin) distmin = dist
if(dist > distmax) distmax = dist
endif
enddo
enddo
enddo
+ print *,'done processing ',NSPEC,' elements out of ',NSPEC
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 *,'minimum length divided in ',NGLLX-1,' average segments (m) = ',distmin/real(NGLLX-1)
print *,'maximum length of an edge in the whole mesh (m) = ',distmax
print *
- end program analyze_CUBIT_Abaqus_mesh
+! ************* 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
+
+ equiangle_skewness_max = - HUGEVAL
+ edge_aspect_ratio_max = - HUGEVAL
+ diagonal_aspect_ratio_max = - HUGEVAL
+ stability_max = - HUGEVAL
+
+! 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)
+ 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)
+ endif
+
+! 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)
+
+ 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)
+
+ 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 *,'edge aspect ratio: 1. perfect above 1. gives stretching factor'
+ print *,'diagonal aspect ratio: 1. perfect above 1. gives stretching factor'
+ print *,'------------'
+
+ print *
+ print *,'max equiangle skewness = ',equiangle_skewness_max
+! print *,'min equiangle skewness = ',equiangle_skewness_min
+ print *
+ print *,'skewness max deviation angle from 90 degrees = ',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
+ print *
+
+! create statistics about mesh quality
+
+ print *
+ print *,'creating histogram and statistics of mesh quality'
+ print *
+
+! 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)
+ 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)
+ 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 *,'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 OpenDX file with elements in a certain range of skewness
+
+ if(USE_OPENDX) then
+
+ print *
+ print *,'creating OpenDX file with subset of elements in skewness range'
+ print *,'between ',skewness_AVS_DX_min,' and ',skewness_AVS_DX_max
+ 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
+ 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)
+ 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)
+ endif
+
+! check if element belongs to requested skewness range
+! and flag all the points to remove multiples
+ if(equiangle_skewness >= skewness_AVS_DX_min .and. equiangle_skewness <= skewness_AVS_DX_max) &
+ ntotspecAVS_DX = ntotspecAVS_DX + 1
+
+ enddo
+
+ if(ntotspecAVS_DX == 0) stop 'no elements in skewness range, no file created'
+
+ 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
+ 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)
+ 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)
+ endif
+
+! check if element belongs to requested skewness range
+! and flag all the points to remove multiples
+ if(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,"(i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6)") &
+ ibool(1,ispec)-1, ibool(4,ispec)-1, ibool(2,ispec)-1, ibool(3,ispec)-1
+ else
+ write(11,"(i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6)") &
+ 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
+ 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 = 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)
+ 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)
+ endif
+
+! check if element belongs to requested skewness range
+! and flag all the points to remove multiples
+ if(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)
+
+ print *
+ print *,'there are ',ntotspecAVS_DX,' elements in AVS or DX skewness range ',skewness_AVS_DX_min,skewness_AVS_DX_max
+ print *
+
+ 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)
+
+ implicit none
+
+ include "constants.h"
+
+ integer iface,icorner,jcorner,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)
+ 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)
+
+ enddo
+ enddo
+
+! compute edge aspect ratio using the NGNOD corners of the element
+ distmin = + HUGEVAL
+ distmax = - HUGEVAL
+ do icorner = 1,NGNOD
+ do jcorner = icorner + 1,NGNOD
+ dist = sqrt((xelm(jcorner) - xelm(icorner))**2 + (yelm(jcorner) - yelm(icorner))**2 + (zelm(jcorner) - zelm(icorner))**2)
+ distmin = min(distmin,dist)
+ distmax = max(distmax,dist)
+ enddo
+ enddo
+ 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)
+ distmin = min(dist1,dist2,dist3,dist4)
+ distmax = max(dist1,dist2,dist3,dist4)
+ diagonal_aspect_ratio = distmax / distmin
+
+ 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)
+
+ implicit none
+
+ include "constants.h"
+
+ integer iface,icorner,jcorner,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(1,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) = 3
+ faces_topo(1,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)
+ equiangle_skewness = - HUGEVAL
+ do iface = 1,1 !! DK DK only one face in 2D
+ 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)
+
+ enddo
+ enddo
+
+! compute edge aspect ratio using the NGNOD corners of the element
+ distmin = + HUGEVAL
+ distmax = - HUGEVAL
+ do icorner = 1,NGNOD
+ do jcorner = icorner + 1,NGNOD
+ dist = sqrt((xelm(jcorner) - xelm(icorner))**2 + (yelm(jcorner) - yelm(icorner))**2 + (zelm(jcorner) - zelm(icorner))**2)
+ distmin = min(distmin,dist)
+ distmax = max(distmax,dist)
+ enddo
+ enddo
+ 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)
+ distmin = min(dist1,dist2)
+ distmax = max(dist1,dist2)
+ diagonal_aspect_ratio = distmax / distmin
+
+ end subroutine create_mesh_quality_data_2D
+
Added: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/constants.h
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/constants.h (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/constants.h 2009-02-24 14:51:55 UTC (rev 14135)
@@ -0,0 +1,12 @@
+
+! number of GLL points in each direction of an element (degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLY = NGLLX
+ integer, parameter :: NGLLZ = NGLLX
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.d+30,TINYVAL = 1.d-9
+
+! some useful constants
+ double precision, parameter :: PI = 3.141592653589793d0
+
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/multiply_CUBIT_Abaqus_mesh_by_1000.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/multiply_CUBIT_Abaqus_mesh_by_1000.f90 2009-02-24 14:50:07 UTC (rev 14134)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/multiply_CUBIT_Abaqus_mesh_by_1000.f90 2009-02-24 14:51:55 UTC (rev 14135)
@@ -8,32 +8,33 @@
implicit none
-! case of HOMOGENE_3D_in_meters.inp
-! integer, parameter :: NPOIN = 21212
+ character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_3D_lisse_300_in_meters.inp'
+ integer, parameter :: NPOIN = 98692
-! case of regolite_3D_rego3d_70m_in_meters.inp
+! character(len=100), parameter :: cubit_mesh_file = 'regolite_3D_rego3d_70m_in_meters.inp'
! integer, parameter :: NPOIN = 4050696
-! case of HOMOGENE_2D_in_meters.inp
+! character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_2D_in_meters.inp'
! integer, parameter :: NPOIN = 3882
-! case of Eros_complet_regolite_fractures_2D_in_meters.inp
- integer, parameter :: NPOIN = 57630
+! character(len=100), parameter :: cubit_mesh_file = 'Eros_complet_regolite_fractures_2D_in_meters.inp'
+! integer, parameter :: NPOIN = 57630
real, dimension(NPOIN) :: x,y,z
integer :: i,iread
! read the mesh
+ open(unit=10,file=cubit_mesh_file,status='unknown',action='read')
! skip the header
- read(*,*)
- read(*,*)
- read(*,*)
+ read(10,*)
+ read(10,*)
+ read(10,*)
! read the points
do i = 1,NPOIN
- read(*,*) iread,x(i),y(i),z(i)
+ 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'
@@ -41,5 +42,7 @@
print *,iread,',',x(i)*1000,',',y(i)*1000,',',z(i)*1000
enddo
+ close(10)
+
end program multiply_CUBIT_Abaqus_mesh
More information about the CIG-COMMITS
mailing list