[cig-commits] r14126 - in seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh: . analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Feb 23 14:23:34 PST 2009
Author: dkomati1
Date: 2009-02-23 14:23:34 -0800 (Mon, 23 Feb 2009)
New Revision: 14126
Added:
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/multiply_CUBIT_Abaqus_mesh_by_1000.f90
Log:
added small programs to analyze a CUBIT mesh stored in ABAQUS ASCII format.
Added: seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/analyze_CUBIT_Abaqus_mesh.f90 2009-02-23 22:23:34 UTC (rev 14126)
@@ -0,0 +1,111 @@
+
+ program analyze_CUBIT_Abaqus_mesh
+
+! 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.
+
+ implicit none
+
+! 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.
+
+! case of regolite_3D_rego3d_70m_in_meters.inp
+! integer, parameter :: NPOIN = 4050696, NELEM = 3410265, NGNOD = 8
+! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+
+! case of HOMOGENE_2D_in_meters.inp
+! integer, parameter :: NPOIN = 3882, NELEM = 3744, NGNOD = 4
+! logical, parameter :: IGNORE_OTHER_HEADERS = .false.
+
+! case of Eros_complet_regolite_fractures_2D_in_meters.inp
+ integer, parameter :: NPOIN = 57630, NELEM = 107904, NGNOD = 4
+ logical, parameter :: IGNORE_OTHER_HEADERS = .true.
+
+ real, dimension(NPOIN) :: x,y,z
+
+ integer, dimension(NGNOD,NELEM) :: ibool
+
+ integer :: i,j,ielem,iread
+
+ real :: distmin,distmax,dist
+
+! read the mesh
+
+ print *,'start reading the CUBIT file'
+
+! skip the header
+ read(*,*)
+ read(*,*)
+ read(*,*)
+
+! read the points
+ do i = 1,NPOIN
+ read(*,*) 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(*,*)
+
+! read the elements
+ do i = 1,NELEM
+
+ 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(*,*)
+
+ read(*,*) 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'
+ endif
+
+ if(iread /= i) then
+ print *,'error at i,iread = ',i,iread
+ stop 'wrong input for an element'
+ endif
+
+ enddo
+
+ 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 ielem = 1,NELEM
+ if(mod(ielem,10000) == 0) print *,'processed ',ielem,' elements out of ',NELEM
+ 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)
+ if(dist < distmin) distmin = dist
+ if(dist > distmax) distmax = dist
+ endif
+ enddo
+ enddo
+ enddo
+
+ 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 *,'maximum length of an edge in the whole mesh (m) = ',distmax
+ print *
+
+ end program analyze_CUBIT_Abaqus_mesh
+
Added: 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 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/multiply_CUBIT_Abaqus_mesh_by_1000.f90 2009-02-23 22:23:34 UTC (rev 14126)
@@ -0,0 +1,45 @@
+
+ program multiply_CUBIT_Abaqus_mesh
+
+! Dimitri Komatitsch, University of Pau, France, February 2009.
+
+! multiply a CUBIT mesh stored in Abaqus ASCII format by 1000
+! to convert it from km to m
+
+ implicit none
+
+! case of HOMOGENE_3D_in_meters.inp
+! integer, parameter :: NPOIN = 21212
+
+! case of regolite_3D_rego3d_70m_in_meters.inp
+! integer, parameter :: NPOIN = 4050696
+
+! case of HOMOGENE_2D_in_meters.inp
+! integer, parameter :: NPOIN = 3882
+
+! case of Eros_complet_regolite_fractures_2D_in_meters.inp
+ integer, parameter :: NPOIN = 57630
+
+ real, dimension(NPOIN) :: x,y,z
+
+ integer :: i,iread
+
+! read the mesh
+
+! skip the header
+ read(*,*)
+ read(*,*)
+ read(*,*)
+
+! read the points
+ do i = 1,NPOIN
+ read(*,*) 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
+ print *,iread,',',x(i)*1000,',',y(i)*1000,',',z(i)*1000
+ enddo
+
+ end program multiply_CUBIT_Abaqus_mesh
+
More information about the CIG-COMMITS
mailing list