[cig-commits] r19361 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES: . tpv16.crack tpv16.crack/description
surendra at geodynamics.org
surendra at geodynamics.org
Sat Jan 14 03:33:25 PST 2012
Author: surendra
Date: 2012-01-14 03:33:25 -0800 (Sat, 14 Jan 2012)
New Revision: 19361
Added:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.jou
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/description/
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/description/TPV16_17_Description_v03.pdf
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/interp5GLL.m
Log:
Added input files for TPV16
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.jou
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.jou (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.jou 2012-01-14 11:33:25 UTC (rev 19361)
@@ -0,0 +1,114 @@
+# Parameters for smaller mesh.
+#
+#{Units('si')}
+#
+# ----------------------------------------------------------------------
+# Parameters for geometry.
+# ----------------------------------------------------------------------
+# Faults have a common width of 15 km
+#{faultWidth=19.5*km}
+#
+# Main fault is 28.0 km long.
+#{mainFaultLength=48.0*km}
+#{mainFaultLength1=48.0*km}
+#
+#
+#
+# ----------------------------------------------------------------------
+# Parameters for mesh.
+# ----------------------------------------------------------------------
+#{dx=4*75*m}
+#
+#{bias_factor=1.1}
+
+
+
+
+
+
+
+# ----------------------------------------------------------------------
+# Geometry
+# ----------------------------------------------------------------------
+reset
+
+create cylinder z {blockHeight} radius {mainFaultLength1}
+volume 1 move z {-blockHeight/2.0}
+
+create planar surface with plane xplane offset 0
+webcut volume 1 with sheet body 2
+delete body 2
+
+move volume 1 x {mainFaultLength/2.0}
+move volume 3 x {-mainFaultLength/2.0}
+
+brick x {mainFaultLength} y {2.0*mainFaultLength1} z {blockHeight}
+volume 4 move z {-blockHeight/2.0}
+
+webcut volume 4 with zplane offset {-faultWidth}
+webcut volume 1 3 with zplane offset {-faultWidth}
+
+webcut volume 4 with yplane
+webcut volume 1 3 with yplane
+
+surface 41 46 merge off
+
+merge tolerance 1e-0
+
+imprint all
+merge all
+
+# ----------------------------------------------------------------------
+# Generate the mesh
+# ----------------------------------------------------------------------
+
+surface 41 46 size {dx}
+mesh surface 41 46
+
+curve 76 scheme bias fine size {dx} factor {bias_factor} start vertex 42
+curve 78 scheme bias fine size {dx} factor {bias_factor} start vertex 41
+mesh surface 50
+curve 77 scheme bias fine size {dx} factor {bias_factor} start vertex 39
+curve 75 scheme bias fine size {dx} factor {bias_factor} start vertex 40
+mesh volume 8
+
+copy mesh volume 8 onto volume 4 nosmoothing
+
+curve 43 scheme bias fine size {dx} factor {bias_factor} start vertex 24
+curve 44 scheme bias fine size {dx} factor {bias_factor} start vertex 26
+curve 45 scheme bias fine size {dx} factor {bias_factor} start vertex 23
+curve 46 scheme bias fine size {dx} factor {bias_factor} start vertex 25
+copy mesh curve 69 77 onto curve 14 source combine
+copy mesh curve 69 77 onto curve 10 source combine
+copy mesh curve 31 onto curve 24
+copy mesh curve 31 onto curve 26
+mesh surface 14
+mesh volume 5
+
+mesh volume 6 7
+
+mesh volume 9 10
+
+copy mesh volume 9 onto volume 1 nosmoothing
+copy mesh volume 9 onto volume 3 nosmoothing
+
+
+########
+
+group "uppr" add node in surface 41
+group "uppr" remove node with x_coord = {-mainFaultLength/2.0}
+group "uppr" remove node with x_coord = {mainFaultLength/2.0}
+group "uppr" remove node with z_coord = {-faultWidth}
+nodeset 400 group uppr
+
+
+group "lowr" add node in surface 46
+group "lowr" remove node with x_coord = {-mainFaultLength/2.0}
+group "lowr" remove node with x_coord = {mainFaultLength/2.0}
+group "lowr" remove node with z_coord = {-faultWidth}
+nodeset 401 group lowr
+
+nodeset move 400 x 0.01
+nodeset move 401 x -0.01
+
+
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.py (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.py 2012-01-14 11:33:25 UTC (rev 19361)
@@ -0,0 +1,147 @@
+#!python
+#!/usr/bin/env python
+
+import cubit
+import cubit2specfem3d
+
+import os
+import sys
+from save_fault_nodes_elements import *
+from absorbing_boundary import *
+
+cubit.cmd('playback "TPV16.jou" ')
+
+os.system('mkdir -p MESH')
+
+
+
+Au = [41] # A_up
+Ad = [46] # A_down
+
+faultA = fault_input(1,Au,Ad)
+quads_Aup,quads_Adp = save_cracks(faultA.name,faultA.surface_u,faultA.surface_d)
+#Unpacking list.
+quads_Au=unpack_list(quads_Aup)
+quads_Ad=unpack_list(quads_Adp)
+
+print 'len(Au):',len(quads_Au)
+print 'len(Ad):',len(quads_Ad)
+
+if not (len(quads_Au)==len(quads_Ad)):
+ print 'Number of elements for each fauld side up and down do not concide'
+
+save_elements_nodes(faultA.name,quads_Au,quads_Ad)
+
+
+xmin = [63,69,40]
+xmax = [34,55,57]
+ymin = [22,26]
+ymax = [20,28]
+zmin = [10,14,8]
+zmax = [62,42,52,60,50,70]
+entities=['face']
+define_bc(entities,xmin,xmax,ymin,ymax,zmin,zmax)
+
+cubit.cmd('block 1 name "elastic 1" ') # material region
+cubit.cmd('block 1 attribute count 5')
+cubit.cmd('block 1 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 1 attribute index 2 6000') # vp
+cubit.cmd('block 1 attribute index 3 3464') # vs
+cubit.cmd('block 1 attribute index 4 2670') # rho
+cubit.cmd('block 1 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+cubit.cmd('block 2 name "elastic 2" ') # material region
+cubit.cmd('block 2 attribute count 5')
+cubit.cmd('block 2 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 2 attribute index 2 6000') # vp
+cubit.cmd('block 2 attribute index 3 3464') # vs
+cubit.cmd('block 2 attribute index 4 2670') # rho
+cubit.cmd('block 2 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+cubit.cmd('block 3 name "elastic 3" ') # material region
+cubit.cmd('block 3 attribute count 5')
+cubit.cmd('block 3 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 3 attribute index 2 6000') # vp
+cubit.cmd('block 3 attribute index 3 3464') # vs
+cubit.cmd('block 3 attribute index 4 2670') # rho
+cubit.cmd('block 3 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+
+cubit.cmd('block 4 name "elastic 4" ') # material region
+cubit.cmd('block 4 attribute count 5')
+cubit.cmd('block 4 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 4 attribute index 2 6000') # vp
+cubit.cmd('block 4 attribute index 3 3464') # vs
+cubit.cmd('block 4 attribute index 4 2670') # rho
+cubit.cmd('block 4 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+
+
+cubit.cmd('block 5 name "elastic 5" ') # material region
+cubit.cmd('block 5 attribute count 5')
+cubit.cmd('block 5 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 5 attribute index 2 6000') # vp
+cubit.cmd('block 5 attribute index 3 3464') # vs
+cubit.cmd('block 5 attribute index 4 2670') # rho
+cubit.cmd('block 5 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+
+
+cubit.cmd('block 6 name "elastic 6" ') # material region
+cubit.cmd('block 6 attribute count 5')
+cubit.cmd('block 6 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 6 attribute index 2 6000') # vp
+cubit.cmd('block 6 attribute index 3 3464') # vs
+cubit.cmd('block 6 attribute index 4 2670') # rho
+cubit.cmd('block 6 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+
+
+cubit.cmd('block 7 name "elastic 7" ') # material region
+cubit.cmd('block 7 attribute count 5')
+cubit.cmd('block 7 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 7 attribute index 2 6000') # vp
+cubit.cmd('block 7 attribute index 3 3464') # vs
+cubit.cmd('block 7 attribute index 4 2670') # rho
+cubit.cmd('block 7 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+
+
+cubit.cmd('block 8 name "elastic 8" ') # material region
+cubit.cmd('block 8 attribute count 5')
+cubit.cmd('block 8 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 8 attribute index 2 6000') # vp
+cubit.cmd('block 8 attribute index 3 3464') # vs
+cubit.cmd('block 8 attribute index 4 2670') # rho
+cubit.cmd('block 8 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+
+
+cubit.cmd('block 9 name "elastic 9" ') # material region
+cubit.cmd('block 9 attribute count 5')
+cubit.cmd('block 9 attribute index 1 1') # flag for fault side 1
+cubit.cmd('block 9 attribute index 2 6000') # vp
+cubit.cmd('block 9 attribute index 3 3464') # vs
+cubit.cmd('block 9 attribute index 4 2670') # rho
+cubit.cmd('block 9 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+
+
+#cubit.cmd('block 10 name "elastic 10" ') # material region
+#cubit.cmd('block 10 attribute count 5')
+#cubit.cmd('block 10 attribute index 1 1') # flag for fault side 1
+#cubit.cmd('block 10 attribute index 2 6000') # vp
+#cubit.cmd('block 10 attribute index 3 3464') # vs
+#cubit.cmd('block 10 attribute index 4 2670') # rho
+#cubit.cmd('block 10 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+
+
+
+#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
+
+cubit2specfem3d.export2SESAME('MESH')
+
+
+
+
+
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/description/TPV16_17_Description_v03.pdf
===================================================================
(Binary files differ)
Property changes on: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/description/TPV16_17_Description_v03.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/interp5GLL.m
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/interp5GLL.m (rev 0)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/interp5GLL.m 2012-01-14 11:33:25 UTC (rev 19361)
@@ -0,0 +1,129 @@
+data = load('tpv17_input_file.txt');
+
+ ix = data(:,1)+1;
+ iy = data(:,2)+1;
+ xx = data(:,3);
+ yy = data(:,4);
+ zz5 = data(:,5);
+ zz6 = data(:,6); % initial shear stress
+ zz7 = data(:,7);
+ zz8 = data(:,8);
+ zz9 = data(:,9);
+ zz10 = data(:,10);
+ zz11 = data(:,11);
+ zz12 = data(:,12);
+ zz13 = data(:,13);
+ zz14 = data(:,14);
+
+ for i=1:length(ix)
+ X(ix(i),iy(i)) = xx(i);
+ Y(ix(i),iy(i)) = yy(i);
+ Z5(ix(i),iy(i)) = zz5(i);
+ Z6(ix(i),iy(i)) = zz6(i);
+ Z7(ix(i),iy(i)) = zz7(i);
+ Z8(ix(i),iy(i)) = zz8(i);
+ Z9(ix(i),iy(i)) = zz9(i);
+ Z10(ix(i),iy(i)) = zz10(i);
+ Z11(ix(i),iy(i)) = zz11(i);
+ Z12(ix(i),iy(i)) = zz12(i);
+ Z13(ix(i),iy(i)) = zz13(i);
+ Z14(ix(i),iy(i)) = zz14(i);
+ end
+
+ figure;
+ subplot(2,1,1);
+ pcolor(X,Y,Z6); shading flat,colorbar,axis image;
+ %orient tall, wysiwyg;
+
+ %**** Set here the parameters of the square box domain and mesh : ****
+ LX = 48e3;
+ LY = 19.5e3;
+ NELX = 160;
+ NELY = 65;
+ P = 4;
+
+ dxe = LX/NELX;
+ dye = LY/NELY;
+ NEL = NELX*NELY;
+ NGLL = P+1; % number of GLL nodes per element
+
+ %[iglob,x,y] = MeshBox(LX,LY,NELX,NELY,NGLL);
+
+ XGLL = GetGLL(NGLL); % local cordinate of GLL quadrature points
+
+ iglob = zeros(NGLL,NGLL,NELX*NELY); % local to global index mapping
+ nglob = (NELX*(NGLL-1)+1)*(NELY*(NGLL-1)+1); % number of global nodes
+
+ x = zeros(nglob,1); % coordinates of GLL nodes
+ y = zeros(nglob,1);
+
+ e=0;
+ last_iglob = 0;
+ igL = reshape([1:NGLL*(NGLL-1)],NGLL-1,NGLL);
+ igB = reshape([1:NGLL*(NGLL-1)],NGLL,NGLL-1);
+ igLB = reshape([1:(NGLL-1)*(NGLL-1)],NGLL-1,NGLL-1);
+ xgll = repmat( 0.5*(1+XGLL) , 1,NGLL);
+ ygll = dye*xgll';
+ xgll = dxe*xgll;
+
+ for ey=1:NELY,
+ for ex=1:NELX,
+ e = e+1;
+ % Take care of redundant nodes at element edges :
+ if e==1 % first element: bottom-left
+ ig = reshape([1:NGLL*NGLL],NGLL,NGLL);
+ else
+ if ey==1 % elements on first (bottom) row
+ ig(1,:) = iglob(NGLL,:,e-1); % left edge
+ ig(2:NGLL,:) = last_iglob + igL; % the rest
+ elseif ex==1 % elements on first (left) column
+ ig(:,1) = iglob(:,NGLL,e-NELX); % bottom edge
+ ig(:,2:NGLL) = last_iglob + igB; % the rest
+ else % other elements
+ ig(1,:) = iglob(NGLL,:,e-1); % left edge
+ ig(:,1) = iglob(:,NGLL,e-NELX); % bottom edge
+ ig(2:NGLL,2:NGLL) = last_iglob + igLB;
+ end
+ end
+ iglob(:,:,e) = ig;
+ last_iglob = ig(NGLL,NGLL);
+
+ % Global coordinates of the computational (GLL) nodes
+ x(ig) = dxe*(ex-1)+xgll;
+ y(ig) = dye*(ey-1)+ygll;
+
+ end
+ end
+
+ [ys I]=sort(y);
+ xs = x(I);
+ for i=1:length(x)
+ XSEM(ix(i),iy(i)) = xs(i);
+ YSEM(ix(i),iy(i)) = ys(i);
+ end
+
+ ZSEM5 = griddata(X,Y,Z5,XSEM,YSEM,'linear');
+ ZSEM6 = griddata(X,Y,Z6,XSEM,YSEM,'linear');
+ ZSEM7 = griddata(X,Y,Z7,XSEM,YSEM,'linear');
+ ZSEM8 = griddata(X,Y,Z8,XSEM,YSEM,'linear');
+ ZSEM9 = griddata(X,Y,Z9,XSEM,YSEM,'linear');
+ ZSEM10 = griddata(X,Y,Z10,XSEM,YSEM,'linear');
+ ZSEM11 = griddata(X,Y,Z11,XSEM,YSEM,'linear');
+ ZSEM12 = griddata(X,Y,Z12,XSEM,YSEM,'linear');
+ ZSEM13 = griddata(X,Y,Z13,XSEM,YSEM,'linear');
+ ZSEM14 = griddata(X,Y,Z14,XSEM,YSEM,'linear');
+
+ subplot(2,1,2);
+ pcolor(XSEM,YSEM,ZSEM6); shading flat,colorbar,axis image;
+
+ file1 = sprintf('tpv17_input_file_interpGLL5.txt');
+ fid = fopen(file1,'w');
+ for i=1:length(ix)
+ fprintf(fid,'%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e%14.6e\n',...
+ ix(i)-1,iy(i)-1,XSEM(ix(i),iy(i)),YSEM(ix(i),iy(i)),ZSEM5(ix(i),iy(i)),ZSEM6(ix(i),iy(i)),...
+ ZSEM7(ix(i),iy(i)),ZSEM8(ix(i),iy(i)),ZSEM9(ix(i),iy(i)),ZSEM10(ix(i),iy(i)),ZSEM11(ix(i),iy(i)),...
+ ZSEM12(ix(i),iy(i)),ZSEM13(ix(i),iy(i)),ZSEM14(ix(i),iy(i)));
+ end
+ fclose(fid);
+
+ dataSEM = load('tpv17_input_file_interpGLL5.txt');
More information about the CIG-COMMITS
mailing list