[cig-commits] r14229 - seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Mar 4 12:59:07 PST 2009
Author: dkomati1
Date: 2009-03-04 12:59:07 -0800 (Wed, 04 Mar 2009)
New Revision: 14229
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90
Log:
fixed an important bug in the calculation of the minimum and maximum edge length
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-04 17:27:56 UTC (rev 14228)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 2009-03-04 20:59:07 UTC (rev 14229)
@@ -26,7 +26,7 @@
! 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.
+! Dimitri Komatitsch, University of Pau, France, March 2009.
!! DK DK
!! DK DK this routine could be improved:
@@ -44,56 +44,58 @@
! 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 = 'mesh_piero_doubling.inp'
-! integer, parameter :: NPOIN = 4044574, NSPEC = 3912880, NGNOD = 8
-! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
-! double precision, parameter :: delta_t = 1.d-3
+ 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
-
! 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 = 0.15e-3
+! 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
- 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,j,ispec,iread,iformat
+ integer :: i,ispec,iread,iformat,ispec_min_edge_length,ispec_max_edge_length
- double precision :: distmin,distmax,dist
+ 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
+ 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
! for stability
- double precision stability
- double precision stability_min
- double precision stability_max
+ double precision :: stability,stability_min,stability_max
! for histogram
integer, parameter :: NCLASS = 20
@@ -187,33 +189,6 @@
print *,'end reading the CUBIT file'
print *,'start computing the minimum and maximum edge size'
-! compute minimum and maximum distance
- distmin = + 1.e30
- distmax = - 1.e30
-
-! loop on all the elements
- do 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,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 ',NGLLX-1,' average segments (m) = ',distmin/real(NGLLX-1)
- print *,'maximum length of an edge in the whole mesh (m) = ',distmax
- print *
-
! ************* compute min and max of skewness and ratios ******************
! erase minimum and maximum of quality numbers
@@ -221,36 +196,51 @@
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)
+ 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)
+ 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)
+ 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)
+ 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 *
@@ -261,6 +251,10 @@
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 *
@@ -277,13 +271,11 @@
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
@@ -293,10 +285,10 @@
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)
+ 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)
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
endif
! store skewness in histogram
@@ -371,10 +363,10 @@
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)
+ 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)
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
endif
! check if element belongs to requested skewness range
@@ -413,10 +405,10 @@
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)
+ 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)
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
endif
! check if element belongs to requested skewness range
@@ -454,10 +446,10 @@
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)
+ 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)
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
endif
! check if element belongs to requested skewness range
@@ -489,13 +481,13 @@
! 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)
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
implicit none
include "constants.h"
- integer iface,icorner,jcorner,ispec,NSPEC,NPOIN,NGNOD,i
+ integer :: iface,icorner,ispec,NSPEC,NPOIN,NGNOD,i
double precision, dimension(NPOIN) :: x,y,z
@@ -592,7 +584,11 @@
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
@@ -616,31 +612,28 @@
! compute equiangle skewness
equiangle_skewness = max(equiangle_skewness,dabs(2.d0 * angle_vectors - PI) / PI)
- enddo
- enddo
+! 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)
-! 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
+! 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)
- distmin = min(dist1,dist2,dist3,dist4)
- distmax = max(dist1,dist2,dist3,dist4)
- diagonal_aspect_ratio = distmax / distmin
+ 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
@@ -651,13 +644,13 @@
! 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)
+ equiangle_skewness,edge_aspect_ratio,diagonal_aspect_ratio,stability,distmin,distmax)
implicit none
include "constants.h"
- integer icorner,jcorner,ispec,NSPEC,NPOIN,NGNOD,i
+ integer :: icorner,ispec,NSPEC,NPOIN,NGNOD,i
double precision, dimension(NPOIN) :: x,y,z
@@ -725,7 +718,11 @@
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
@@ -748,28 +745,25 @@
! compute equiangle skewness
equiangle_skewness = max(equiangle_skewness,dabs(2.d0 * angle_vectors - PI) / PI)
- enddo
+! 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)
-! 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
+ distmin = min(distmin,dist)
+ distmax = max(distmax,dist)
+
enddo
- edge_aspect_ratio = distmax / distmin
+! 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)
- distmin = min(dist1,dist2)
- distmax = max(dist1,dist2)
- diagonal_aspect_ratio = distmax / distmin
+ 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