[cig-commits] r14250 - seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Mar 7 17:27:34 PST 2009
Author: dkomati1
Date: 2009-03-07 17:27:34 -0800 (Sat, 07 Mar 2009)
New Revision: 14250
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90
Log:
fixed some bugs and added an option to extract a single element
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-07 23:45:33 UTC (rev 14249)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 2009-03-08 01:27:34 UTC (rev 14250)
@@ -24,16 +24,13 @@
!=====================================================================
! 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
+! 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:
+!! DK DK this routine could be improved by computing the mean in addition to min and max of ratios
!! 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
@@ -44,11 +41,11 @@
! 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 = '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
@@ -74,25 +71,25 @@
! 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
+ 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
+ integer :: i,ispec,iread,iformat,ispec_min_edge_length,ispec_max_edge_length, &
+ ispec_begin,ispec_end,ispec_to_output
- double precision :: distmin,distmax
-
! 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
@@ -100,8 +97,8 @@
! for histogram
integer, parameter :: NCLASS = 20
integer classes_skewness(0:NCLASS-1)
- integer iclass
- double precision current_percent,total_percent
+ integer :: iclass
+ double precision :: current_percent,total_percent
! to export elements that have a certain skewness range to OpenDX
integer :: ntotspecAVS_DX
@@ -110,13 +107,16 @@
if(NGNOD /= 4 .and. NGNOD /= 8) stop 'NGNOD should be 4 or 8'
print *
- print *,'1 = output bad elements above a certain skewness threshold in OpenDX format'
- print *,'2 = do not output bad elements'
+ 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>2) stop 'exiting...'
- if(iformat == 1) then
+
+ if(iformat < 1 .or. iformat > 3) stop 'exiting...'
+
+ if(iformat == 1 .or. iformat == 2) then
USE_OPENDX = .true.
else
USE_OPENDX = .false.
@@ -124,6 +124,8 @@
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
@@ -138,9 +140,16 @@
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')
@@ -186,7 +195,9 @@
enddo
close(10)
- print *,'end reading the CUBIT file'
+ print *,'done reading the CUBIT file'
+ print *
+
print *,'start computing the minimum and maximum edge size'
! ************* compute min and max of skewness and ratios ******************
@@ -330,27 +341,35 @@
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
+ if(equiangle_skewness_max >= 0.75d0) then
print *
print *,'*********************************************'
print *,'*********************************************'
- print *,' WARNING, mesh is bad (max skewness >= 0.80)'
+ 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 *
- print *,'creating OpenDX file with subset of elements in skewness range'
- print *,'between ',skewness_AVS_DX_min,' and ',skewness_AVS_DX_max
+ 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 *************
@@ -359,6 +378,8 @@
ntotspecAVS_DX = 0
! loop on all the elements
+ if(iformat == 1) then
+
do ispec = 1,NSPEC
if(NGNOD == 4) then
@@ -370,15 +391,19 @@
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
+ else
+! outputing a single element
+ ntotspecAVS_DX = 1
+ endif
+
if(ntotspecAVS_DX == 0) then
stop 'no elements in skewness range, no file created'
- else
+ 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 *
@@ -401,8 +426,16 @@
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(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)
@@ -411,21 +444,21 @@
equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
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
+! 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,"(i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6)") &
+ 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,"(i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6,1x,i6)") &
+ 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
- print *,'element ',ispec,' belongs to the range and has skewness = ',sngl(equiangle_skewness)
+ if(iformat == 1) print *,'element ',ispec,' belongs to the range and has skewness = ',sngl(equiangle_skewness)
endif
enddo
@@ -442,7 +475,7 @@
write(11,*) 'object 3 class array type float rank 0 items ',ntotspecAVS_DX,' data follows'
! loop on all the elements
- do ispec = 1,NSPEC
+ 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, &
@@ -452,9 +485,9 @@
equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
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) &
+! 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
@@ -613,9 +646,7 @@
equiangle_skewness = max(equiangle_skewness,dabs(2.d0 * angle_vectors - PI) / PI)
! compute min and max size of an edge
- dist = sqrt((xelm(faces_topo(iface,icorner+1)) - xelm(faces_topo(iface,icorner)))**2 + &
- (yelm(faces_topo(iface,icorner+1)) - yelm(faces_topo(iface,icorner)))**2 + &
- (zelm(faces_topo(iface,icorner+1)) - zelm(faces_topo(iface,icorner)))**2)
+ dist = sqrt(vectorA_x**2 + vectorA_y**2 + vectorA_z**2)
distmin = min(distmin,dist)
distmax = max(distmax,dist)
@@ -746,9 +777,7 @@
equiangle_skewness = max(equiangle_skewness,dabs(2.d0 * angle_vectors - PI) / PI)
! compute min and max size of an edge
- dist = sqrt((xelm(faces_topo(icorner+1)) - xelm(faces_topo(icorner)))**2 + &
- (yelm(faces_topo(icorner+1)) - yelm(faces_topo(icorner)))**2 + &
- (zelm(faces_topo(icorner+1)) - zelm(faces_topo(icorner)))**2)
+ dist = sqrt(vectorA_x**2 + vectorA_y**2 + vectorA_z**2)
distmin = min(distmin,dist)
distmax = max(distmax,dist)
More information about the CIG-COMMITS
mailing list