[cig-commits] [commit] devel: First test HEX27 earth chunk mesher (e58732b)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 20 10:08:06 PST 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/b43c63c4cea034e3e6f647cebc591fb4c2ab3894...fb5ff57e01f82037a9b22ebd3c3c184e9e4f14f6
>---------------------------------------------------------------
commit e58732b79c970425e123ea14f700700255116036
Author: Clément Durochat <c.durochat at gmail.com>
Date: Mon Nov 3 15:28:20 2014 +0100
First test HEX27 earth chunk mesher
>---------------------------------------------------------------
e58732b79c970425e123ea14f700700255116036
...EX8_Mesher.f90 => earth_chunk_HEX27_Mesher.f90} | 233 +++++++++++++++------
src/meshfem3D/meshfem3D.f90 | 6 +-
src/meshfem3D/rules.mk | 1 +
3 files changed, 170 insertions(+), 70 deletions(-)
diff --git a/src/meshfem3D/earth_chunk_HEX8_Mesher.f90 b/src/meshfem3D/earth_chunk_HEX27_Mesher.f90
similarity index 86%
copy from src/meshfem3D/earth_chunk_HEX8_Mesher.f90
copy to src/meshfem3D/earth_chunk_HEX27_Mesher.f90
index 194fea7..2638757 100644
--- a/src/meshfem3D/earth_chunk_HEX8_Mesher.f90
+++ b/src/meshfem3D/earth_chunk_HEX27_Mesher.f90
@@ -25,7 +25,7 @@
!
!=====================================================================
- subroutine earth_chunk_HEX8_Mesher(NGNOD)
+ subroutine earth_chunk_HEX27_Mesher(NGNOD)
use constants, only: NGLLX, NGLLY, NGLLZ, NDIM
@@ -34,10 +34,9 @@
!==============================================================================================!
! !
! Singular option of meshfem3D : MESH OF A GLOBE EARTH CHUNK FOR THE INTERFACE DSM-SPECFEM3D !
-! Case of 8 nodes per element (HEX8) !
+! Case of 27 nodes per element (HEX27) !
! !
-! VM, February 2013 !
-! Integrated in meshfem3d by CD, September 2014 !
+! Integrated in meshfem3d by CD, October 2014 !
! !
! WARNING : A local convention is used for the mapping of !
! the cubic sphere (to complete) !
@@ -119,7 +118,7 @@
TINYVAL = 1.d-9
ZERO = 0.d0
- open(49, file=trim(MESH)//'output_mesher_chunk_HEX8.txt')
+ open(49, file=trim(MESH)//'output_mesher_chunk_HEX27.txt')
if (RUN_BENCHMARK) then
!
@@ -196,38 +195,126 @@
iboun(:,:) = .false.
+
+!! MODIF HEX27 LA ==> cf call hex_nodes-----------------------------
+
+ ! corner nodes
+
iaddx(1) = 0
iaddy(1) = 0
iaddz(1) = 0
- iaddx(2) = 1
+ iaddx(2) = 2
iaddy(2) = 0
iaddz(2) = 0
- iaddx(3) = 1
- iaddy(3) = 1
+ iaddx(3) = 2
+ iaddy(3) = 2
iaddz(3) = 0
iaddx(4) = 0
- iaddy(4) = 1
+ iaddy(4) = 2
iaddz(4) = 0
iaddx(5) = 0
iaddy(5) = 0
- iaddz(5) = 1
+ iaddz(5) = 2
- iaddx(6) = 1
+ iaddx(6) = 2
iaddy(6) = 0
- iaddz(6) = 1
+ iaddz(6) = 2
- iaddx(7) = 1
- iaddy(7) = 1
- iaddz(7) = 1
+ iaddx(7) = 2
+ iaddy(7) = 2
+ iaddz(7) = 2
iaddx(8) = 0
- iaddy(8) = 1
- iaddz(8) = 1
+ iaddy(8) = 2
+ iaddz(8) = 2
+
+! midside nodes (nodes located in the middle of an edge)
+
+ iaddx(9) = 1
+ iaddy(9) = 0
+ iaddz(9) = 0
+
+ iaddx(10) = 2
+ iaddy(10) = 1
+ iaddz(10) = 0
+
+ iaddx(11) = 1
+ iaddy(11) = 2
+ iaddz(11) = 0
+
+ iaddx(12) = 0
+ iaddy(12) = 1
+ iaddz(12) = 0
+
+ iaddx(13) = 0
+ iaddy(13) = 0
+ iaddz(13) = 1
+
+ iaddx(14) = 2
+ iaddy(14) = 0
+ iaddz(14) = 1
+
+ iaddx(15) = 2
+ iaddy(15) = 2
+ iaddz(15) = 1
+
+ iaddx(16) = 0
+ iaddy(16) = 2
+ iaddz(16) = 1
+
+ iaddx(17) = 1
+ iaddy(17) = 0
+ iaddz(17) = 2
+
+ iaddx(18) = 2
+ iaddy(18) = 1
+ iaddz(18) = 2
+
+ iaddx(19) = 1
+ iaddy(19) = 2
+ iaddz(19) = 2
+
+ iaddx(20) = 0
+ iaddy(20) = 1
+ iaddz(20) = 2
+! side center nodes (nodes located in the middle of a face)
+
+ iaddx(21) = 1
+ iaddy(21) = 1
+ iaddz(21) = 0
+
+ iaddx(22) = 1
+ iaddy(22) = 0
+ iaddz(22) = 1
+
+ iaddx(23) = 2
+ iaddy(23) = 1
+ iaddz(23) = 1
+
+ iaddx(24) = 1
+ iaddy(24) = 2
+ iaddz(24) = 1
+
+ iaddx(25) = 0
+ iaddy(25) = 1
+ iaddz(25) = 1
+
+ iaddx(26) = 1
+ iaddy(26) = 1
+ iaddz(26) = 2
+
+! center node (barycenter of the eight corners)
+
+ iaddx(27) = 1
+ iaddy(27) = 1
+ iaddz(27) = 1
+
+!! --------------------------------------------
!
!===========================================================================
!
@@ -400,6 +487,8 @@
! 8 vertices of the element ispec
do ia=1,NGNOD
+!! MODIF HEX27 LA -----------------------------
+
i=iaddx(ia)
j=iaddy(ia)
k=iaddz(ia)
@@ -416,13 +505,6 @@
y = 2.d0*ratio_eta-1.d0
y = tan((ANGULAR_WIDTH_ETA_RAD/2.d0) * y)
- !if (ilat==0.and.iz==0) write(49,*) ia,i,ratio_xi
- !mapping cubic sphere (k=5) (Chevrot et al 2012)
- ! (opposite sign)
- !pz = z/dsqrt(1.d0 + y*y + x*x)
- !px = x*pz
- !py = y*pz
- ! mapping which makes a chunk at the North Pole
! mapping cubic sphere (k=6, Chevrot at al 2012, avec -z)
pz= z/dsqrt(1.d0 + y*y + x*x) !(=r/s)
@@ -434,15 +516,10 @@
ygrid(i+1,j+1,k+1,ispec) = py !py
zgrid(i+1,j+1,k+1,ispec) = pz
- !xgrid(i+1,j+1,k+1,ispec) = py ! long
- !ygrid(i+1,j+1,k+1,ispec) = pz ! lat
- !zgrid(i+1,j+1,k+1,ispec) = px ! prof
-
xelm(ia)=xgrid(i+1,j+1,k+1,ispec)
yelm(ia)=ygrid(i+1,j+1,k+1,ispec)
zelm(ia)=zgrid(i+1,j+1,k+1,ispec)
-
enddo
! INTERFACE FOR DSM ------
@@ -507,19 +584,9 @@
if (iz==0) then ! pas besoin du test comme précédemment car je stocke tout dans des tableaux et c'est pas
! grave si on récrit les memes choses
call calc_gll_points(xelm,yelm,zelm,xstore,ystore,zstore,shape3D,NGNOD,NGLLX,NGLLY,NGLLZ)
- call write_Igm_file(42,ispec2Dzmin,NGLLX,NGLLY,ilon,ilat,0,ilayer_current)
- !open(125,file='ggl_elemts')
-!!$ do kkk=1,1!NGLLZ
-!!$ do jjj=1,NGLLY
-!!$ do iii=1,NGLLX
-!!$ write(125,'(3f20.10)') xstore(iii,jjj,kkk)/1000.d0, &
-!!$ ystore(iii,jjj,kkk)/1000.d0, zstore(iii,jjj,kkk)/1000.d0
-!!$ !write(*,*) xstore
-!!$ enddo
-!!$ enddo
-!!$ enddo
- !close(125)
- ! read(*,*) ia
+
+ call write_Igm_file(42,ispec2Dzmin,NGLLX,NGLLY,ilon,ilat,0,ilayer_current)
+
call store_zmin_points(xstore,ystore,zstore,NGLLX,NGLLY,NGLLZ,rotation_matrix,&
lon_zmin,lat_zmin,nlon_dsm,nlat_dsm,ilon,ilat)
endif
@@ -565,9 +632,9 @@
ieoff = 8 * (ispec - 1)
ilocnum = 0
- do k=1,2
- do j=1,2
- do i=1,2
+ do k=1,3
+ do j=1,3
+ do i=1,3
ilocnum = ilocnum + 1
xp(ilocnum + ieoff)= xgrid(i,j,k,ispec)
@@ -585,15 +652,17 @@
deallocate(xp,yp,zp)
allocate(xp(nglob),yp(nglob),zp(nglob))
+!! MODIF HEX27 LA -----------------------------
+
! on ne stocke que les points de la grille et leur numeros
do ispec=1,nspec
ieoff = 8 * (ispec - 1)
ilocnum = 0
- do k=1,2
- do j=1,2
- do i=1,2
- ilocnum=ilocnum+1
- inum_loc(i,j,k,ispec) = iglob(ilocnum+ieoff)
+ do k=1,3
+ do j=1,3
+ do i=1,3
+ ilocnum = ilocnum + 1
+ inum_loc(i,j,k,ispec) = iglob(ilocnum+ieoff)
xp(iglob(ilocnum+ieoff)) = xgrid(i,j,k,ispec)
yp(iglob(ilocnum+ieoff)) = ygrid(i,j,k,ispec)
zp(iglob(ilocnum+ieoff)) = zgrid(i,j,k,ispec)
@@ -634,6 +703,8 @@
!stop
! -------------------------------- SAUVEGARDE DES MESH FILES -----------
+!! MODIF HEX27 LA -----------------------------
+
open(27,file=trim(MESH)//'nodes_coords_file')
write(27,*) nglob ! nb de sommets
do kglob=1,nglob
@@ -644,10 +715,14 @@
open(27,file=trim(MESH)//'mesh_file')
write(27,*) nspec
do ispec=1,nspec
- write(27,'(9i15)') ispec,inum_loc(1,1,1,ispec),inum_loc(2,1,1,ispec),&
- inum_loc(2,2,1,ispec),inum_loc(1,2,1,ispec),&
- inum_loc(1,1,2,ispec),inum_loc(2,1,2,ispec),&
- inum_loc(2,2,2,ispec),inum_loc(1,2,2,ispec)
+ write(27,'(28i15)') ispec, &
+ inum_loc(1,1,1,ispec), inum_loc(3,1,1,ispec), inum_loc(3,3,1,ispec), inum_loc(1,3,1,ispec), &
+ inum_loc(1,1,3,ispec), inum_loc(3,1,3,ispec), inum_loc(3,3,3,ispec), inum_loc(1,3,3,ispec), &
+ inum_loc(2,1,1,ispec), inum_loc(3,2,1,ispec), inum_loc(2,3,1,ispec), inum_loc(1,2,1,ispec), &
+ inum_loc(1,1,2,ispec), inum_loc(3,1,2,ispec), inum_loc(3,3,2,ispec), inum_loc(1,3,2,ispec), &
+ inum_loc(2,1,3,ispec), inum_loc(3,2,3,ispec), inum_loc(2,3,3,ispec), inum_loc(1,2,3,ispec), &
+ inum_loc(2,2,1,ispec), inum_loc(2,1,2,ispec), inum_loc(3,2,2,ispec), inum_loc(2,3,2,ispec), &
+ inum_loc(1,2,2,ispec), inum_loc(2,2,3,ispec), inum_loc(2,2,2,ispec)
enddo
close(27)
!
@@ -655,48 +730,72 @@
open(27,file=trim(MESH)//'absorbing_surface_file_xmin')
write(27,*) ispec2Dxmin
do ispec=1,nspec
- if (iboun(1,ispec)) write(27,'(5(i10,1x))') ispec,inum_loc(1,1,1,ispec),inum_loc(1,2,1,ispec),&
- inum_loc(1,2,2,ispec),inum_loc(1,1,2,ispec)
+ if (iboun(1,ispec)) write(27,'(10(i10,1x))') ispec, &
+ inum_loc(1,1,1,ispec), inum_loc(1,3,1,ispec), &
+ inum_loc(1,3,3,ispec), inum_loc(1,1,3,ispec), &
+ inum_loc(1,2,1,ispec), inum_loc(1,3,2,ispec), &
+ inum_loc(1,2,3,ispec), inum_loc(1,1,2,ispec), &
+ inum_loc(1,2,2,ispec)
enddo
close(27)
open(27,file=trim(MESH)//'absorbing_surface_file_xmax')
write(27,*) ispec2Dxmax
do ispec=1,nspec
- if (iboun(2,ispec)) write(27,'(5(i10,1x))') ispec,inum_loc(2,1,1,ispec),inum_loc(2,2,1,ispec),&
- inum_loc(2,2,2,ispec),inum_loc(2,1,2,ispec)
+ if (iboun(2,ispec)) write(27,'(10(i10,1x))') ispec, &
+ inum_loc(3,1,1,ispec), inum_loc(3,3,1,ispec), &
+ inum_loc(3,3,3,ispec), inum_loc(3,1,3,ispec), &
+ inum_loc(3,2,1,ispec), inum_loc(3,3,2,ispec), &
+ inum_loc(3,2,3,ispec), inum_loc(3,1,2,ispec), &
+ inum_loc(3,2,2,ispec)
enddo
close(27)
open(27,file=trim(MESH)//'absorbing_surface_file_ymin')
write(27,*) ispec2Dymin
do ispec=1,nspec
- if (iboun(3,ispec)) write(27,'(5(i10,1x))') ispec,inum_loc(1,1,1,ispec),inum_loc(2,1,1,ispec),&
- inum_loc(2,1,2,ispec),inum_loc(1,1,2,ispec)
+ if (iboun(3,ispec)) write(27,'(10(i10,1x))') ispec, &
+ inum_loc(1,1,1,ispec), inum_loc(3,1,1,ispec), &
+ inum_loc(3,1,3,ispec), inum_loc(1,1,3,ispec), &
+ inum_loc(2,1,1,ispec), inum_loc(3,1,2,ispec), &
+ inum_loc(2,1,3,ispec), inum_loc(1,1,2,ispec), &
+ inum_loc(2,1,2,ispec)
enddo
close(27)
open(27,file=trim(MESH)//'absorbing_surface_file_ymax')
write(27,*) ispec2Dymax
do ispec=1,nspec
- if (iboun(4,ispec)) write(27,'(5(i10,1x))') ispec,inum_loc(1,2,1,ispec),inum_loc(2,2,1,ispec),&
- inum_loc(2,2,2,ispec),inum_loc(1,2,2,ispec)
+ if (iboun(4,ispec)) write(27,'(10(i10,1x))') ispec, &
+ inum_loc(1,3,1,ispec), inum_loc(3,3,1,ispec), &
+ inum_loc(3,3,3,ispec), inum_loc(1,3,3,ispec), &
+ inum_loc(2,3,1,ispec), inum_loc(3,3,2,ispec), &
+ inum_loc(2,3,3,ispec), inum_loc(1,3,2,ispec), &
+ inum_loc(2,3,2,ispec)
enddo
close(27)
open(27,file=trim(MESH)//'absorbing_surface_file_bottom')
write(27,*) ispec2Dzmin
do ispec=1,nspec
- if (iboun(5,ispec)) write(27,'(5(i10,1x))') ispec,inum_loc(1,1,1,ispec),inum_loc(1,2,1,ispec),&
- inum_loc(2,2,1,ispec),inum_loc(2,1,1,ispec)
+ if (iboun(5,ispec)) write(27,'(10(i10,1x))') ispec, &
+ inum_loc(1,1,1,ispec), inum_loc(3,1,1,ispec), &
+ inum_loc(3,3,1,ispec), inum_loc(1,3,1,ispec), &
+ inum_loc(2,1,1,ispec), inum_loc(3,2,1,ispec), &
+ inum_loc(2,3,1,ispec), inum_loc(1,2,1,ispec), &
+ inum_loc(2,2,1,ispec)
enddo
close(27)
open(27,file=trim(MESH)//'free_surface')
write(27,*) ispec2Dzmax
do ispec=1,nspec
- if (iboun(1,ispec)) write(27,'(5(i10,1x))') ispec,inum_loc(1,1,2,ispec),inum_loc(1,2,2,ispec),&
- inum_loc(2,2,2,ispec),inum_loc(2,1,2,ispec)
+ if (iboun(1,ispec)) write(27,'(10(i10,1x))') ispec, &
+ inum_loc(1,1,3,ispec), inum_loc(3,1,3,ispec), &
+ inum_loc(3,3,3,ispec), inum_loc(1,3,3,ispec), &
+ inum_loc(2,1,3,ispec), inum_loc(3,2,3,ispec), &
+ inum_loc(2,3,3,ispec), inum_loc(1,2,3,ispec), &
+ inum_loc(2,2,3,ispec)
enddo
close(27)
@@ -708,7 +807,7 @@
! nothing to do anymore... bailing out
stop
- end subroutine earth_chunk_HEX8_Mesher
+ end subroutine earth_chunk_HEX27_Mesher
!=======================================================================================================
!
diff --git a/src/meshfem3D/meshfem3D.f90 b/src/meshfem3D/meshfem3D.f90
index b7c8c28..5ef3825 100644
--- a/src/meshfem3D/meshfem3D.f90
+++ b/src/meshfem3D/meshfem3D.f90
@@ -420,10 +420,10 @@
else if (NGNOD == 27) then
! creates mesh in MESH/
-! to uncomment when finished ==> ! call earth_chunk_HEX27_Mesher(NGNOD)
+ call earth_chunk_HEX27_Mesher(NGNOD)
! done with mesher
- stop 'Creating a chunk of the earth Mesh not implmented yet with HEX27 elements'
-! to uncomment when finished ==> ! stop 'Done creating a chunk of the earth Mesh (HEX27 elements), see directory MESH/'
+!! stop 'Creating a chunk of the earth Mesh not implmented yet with HEX27 elements'
+ stop 'Done creating a chunk of the earth Mesh (HEX27 elements), see directory MESH/'
else
stop 'Bad number of nodes per hexahedron : NGNOD must be equal to 8 or 27'
diff --git a/src/meshfem3D/rules.mk b/src/meshfem3D/rules.mk
index 7a42a3c..4fc620f 100644
--- a/src/meshfem3D/rules.mk
+++ b/src/meshfem3D/rules.mk
@@ -43,6 +43,7 @@ meshfem3D_TARGETS = \
meshfem3D_OBJECTS = \
$O/check_mesh_quality.mesh.o \
$O/earth_chunk_HEX8_Mesher.mesh.o \
+ $O/earth_chunk_HEX27_Mesher.mesh.o \
$O/earth_chunk_ReadIasp91.mesh.o \
$O/compute_parameters.mesh.o \
$O/create_regions_mesh.mesh.o \
More information about the CIG-COMMITS
mailing list