[cig-commits] r14139 - seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Feb 24 15:51:23 PST 2009
Author: dkomati1
Date: 2009-02-24 15:51:23 -0800 (Tue, 24 Feb 2009)
New Revision: 14139
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 support for Piero Basini's Antarctica CUBIT model
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 23:45:33 UTC (rev 14138)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/check_mesh_quality_CUBIT_Abaqus.f90 2009-02-24 23:51:23 UTC (rev 14139)
@@ -44,6 +44,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 = '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.
@@ -643,7 +648,7 @@
include "constants.h"
- integer iface,icorner,jcorner,ispec,NSPEC,NPOIN,NGNOD,i
+ integer icorner,jcorner,ispec,NSPEC,NPOIN,NGNOD,i
double precision, dimension(NPOIN) :: x,y,z
@@ -666,7 +671,7 @@
! topology of faces of cube for skewness
! only one face in 2D
- integer faces_topo(1,6)
+ integer faces_topo(6)
! store the corners of this element for the skewness routine
do i = 1,NGNOD
@@ -700,42 +705,40 @@
! 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
+! 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)
+ 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
+ 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))
+ 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(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))
+ 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)
+ 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))
+ 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)
+ 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
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 23:45:33 UTC (rev 14138)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/UTILS/external_mesh/analyze_CUBIT_Abaqus_mesh/multiply_CUBIT_Abaqus_mesh_by_1000.f90 2009-02-24 23:51:23 UTC (rev 14139)
@@ -8,8 +8,8 @@
implicit none
- character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_3D_lisse_300_in_meters.inp'
- integer, parameter :: NPOIN = 98692
+! character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_3D_lisse_300_in_meters.inp'
+! integer, parameter :: NPOIN = 98692
! character(len=100), parameter :: cubit_mesh_file = 'regolite_3D_rego3d_70m_in_meters.inp'
! integer, parameter :: NPOIN = 4050696
@@ -17,8 +17,8 @@
! character(len=100), parameter :: cubit_mesh_file = 'HOMOGENE_2D_in_meters.inp'
! integer, parameter :: NPOIN = 3882
-! character(len=100), parameter :: cubit_mesh_file = 'Eros_complet_regolite_fractures_2D_in_meters.inp'
-! integer, parameter :: NPOIN = 57630
+ character(len=100), parameter :: cubit_mesh_file = 'eros_complexe_2d_regolite_fractures_in_meters.inp'
+ integer, parameter :: NPOIN = 57807
real, dimension(NPOIN) :: x,y,z
More information about the CIG-COMMITS
mailing list