[cig-commits] r4013 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed Jul 12 14:51:22 PDT 2006
Author: willic3
Date: 2006-07-12 14:51:22 -0700 (Wed, 12 Jul 2006)
New Revision: 4013
Added:
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/README.utils
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/faultcalc.par
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readnetgen.par
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.f
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.par
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/scaleucd.f
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/makefile
Log:
Added utilities and sample parameter files, as well as a README
containing brief code descriptions.
The scaling codes are simple and appear to work OK.
The readucd2 code has been tested on a single example (BM5 created with
LaGriT) and appears to work OK, but more testing is probably needed.
I still need to make this whole directory part of the GNU build
procedure and install the executables.
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/README.utils
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/README.utils 2006-07-12 20:36:01 UTC (rev 4012)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/README.utils 2006-07-12 21:51:22 UTC (rev 4013)
@@ -0,0 +1,46 @@
+This directory contains a number of simple fortran utility codes.
+The codes are as follows:
+
+addflt.f: Code to produce nodal coordinates with 'fake' nodes
+ inserted along faults. Requires UCD file + split or
+ slippery node files. This code is not yet finished or
+ tested.
+faultcalc.f: Code to take output from readnetgen, readucd, or
+ readucd2 and create a fault slip distribution
+ according to a set of given polynomial coefficients.
+ This code requires a parameter file as well as a nodal
+ coordinate file and the initial split node input file
+ created by one of the codes listed above (*.fbc file).
+faultcalc2.f: Code identical to faultcalc.f code, except that this
+ version uses the *.fcoord file rather than the nodal
+ coordinate file to get fault coordinate information and
+ uses essentially no memory compared to faultcalc.
+fixqt.f: Simple code to insure that all non-vertex nodes for a
+ quadratic hex mesh occur at midside. Requires
+ coordinate and connectivity files with no comments or
+ units.
+lh2qh.f: Simple code to convert a linear hex mesh to a quadratic
+ (20-node) hex mesh. Requires coordinate and
+ connectivity files with no comments or units.
+lt2qt.f: Simple code to convert a linear tet mesh to a quadratic
+ (10-node) tet mesh. Requires coordinate and
+ connectivity files with no comments or units.
+pylith2ucd.f: Simple code to convert PyLith/LithoMop coordinate and
+ connectivity files into AVS UCD format.
+readnetgen.f: Simple code to translate NetGen neutral format to a
+ format suitable for PyLith/LithoMop. The code requires
+ a parameter file in addition to the NetGen file and
+ produces coordinates, connectivities, boundary
+ conditions, and split/slippery node input files.
+readucd.f: Simple code to translate UCD files from LaGriT to a
+ format suitable for PyLith/LithoMop. The code requires
+ a parameter file in addition to the UCD file and
+ produces coordinates, connectivities, boundary
+ conditions, and split/slippery node input files.
+readucd2.f: Similar to readucd above but with more features and
+ flexibility.
+scalecoord.f: Simple code to scale NetGen coordinates.
+scaleucd.f: Simple code to scale UCD coordinates.
+tetcmp.f: Debugging tool to compute determinant for a tet element.
+ucd3d.f: Obsolete code that was used to translate old TECTON
+ output to UCD format.
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/faultcalc.par
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/faultcalc.par 2006-07-12 20:36:01 UTC (rev 4012)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/faultcalc.par 2006-07-12 21:51:22 UTC (rev 4013)
@@ -0,0 +1,92 @@
+#
+# Parameter file for use with program faultcalc.
+#
+# **********************************************************************
+#
+# Global information:
+#
+# nsurf = the number of surfaces for which to compute split node
+# displacements
+# eps = epsilon value to determine whether displacements are zero.
+# If so, the node is not included as a split node.
+#
+# nsurf eps
+#-----------------------------------------------------------------------
+ 4 1.0d-10
+#-----------------------------------------------------------------------
+#
+# For each surface, enter the following information:
+#
+# smin(j) = minimum (x,y,z) coordinate values defining the surface
+# smax(j) = maximum (x,y,z) coordinate values defining the surface
+# scoef(j,k) = coefficient values defining slip for each direction
+# These coefficients define a general equation of
+# second degree. The slip for each coordinate
+# direction is computed as:
+#
+# s(x,y,z) = s1*x^2 + 2*s2*x*y + s3*y^2 + 2*s4*x*z
+# + 2*s5*y*z + s6*z^2 + 2*s7*x + 2*s8*y
+# + 2*s9*z +s10
+#
+#-----------------------------------------------------------------------
+ 4.0d0 16.0d0
+ 0.0d0 12.0d0
+ -16.0d0 0.0d0
+# (x,y,z) extents for fault portion with uniform slip
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.3535d0 0.0d0 0.3535d0
+# Coefficients for uniform slip
+ 4.0d0 16.0d0
+ 12.0d0 16.0d0
+ -16.0d0 0.0d0
+# (x,y,z) extents for fault portion with linear y-taper
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ -0.0441875d0 0.0d0 -0.0441875d0
+ 0.0d0 0.0d0 0.0d0
+ 1.414d0 0.0d0 1.414d0
+# Coefficients for fault portion with linear y-taper
+ 16.0d0 20.0d0
+ 0.0d0 12.0d0
+ -16.0d0 0.0d0
+# (x,y,z) extents for fault portion with linear x-taper
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ -0.0441875d0 0.0d0 -0.0441875d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 1.7675d0 0.0d0 1.7675d0
+# Coefficients for fault portion with linear x-taper
+ 16.0d0 20.0d0
+ 12.0d0 16.0d0
+ -16.0d0 0.0d0
+# (x,y,z) extents for fault portion with bilinear xy-taper
+ 0.0d0 0.0d0 0.0d0
+ 0.011046875d0 0.0d0 0.011046875d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ 0.0d0 0.0d0 0.0d0
+ -0.17675d0 0.0d0 -0.17675d0
+ -0.2209375d0 0.0d0 -0.2209375d0
+ 0.0d0 0.0d0 0.0d0
+ 7.07d0 0.0d0 7.07d0
+# Coefficients for fault portion with bilinear xy-taper
+#-----------------------------------------------------------------------
Property changes on: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/faultcalc.par
___________________________________________________________________
Name: svn:executable
+ *
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/makefile
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/makefile 2006-07-12 20:36:01 UTC (rev 4012)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/makefile 2006-07-12 21:51:22 UTC (rev 4013)
@@ -65,8 +65,23 @@
readucd: readucd.o
$(FCOMPL) $(opt) -o readucd readucd.o
+readucd2: readucd2.o
+ $(FCOMPL) $(opt) -o readucd2 readucd2.o
+
tetcmp: tetcmp.o
$(FCOMPL) $(opt) -o tetcmp tetcmp.o
+faultcalc: faultcalc.o
+ $(FCOMPL) $(opt) -o faultcalc faultcalc.o
+
+faultcalc2: faultcalc2.o
+ $(FCOMPL) $(opt) -o faultcalc2 faultcalc2.o
+
+scalecoord: scalecoord.o
+ $(FCOMPL) $(opt) -o scalecoord scalecoord.o
+
+scaleucd: scaleucd.o
+ $(FCOMPL) $(opt) -o scaleucd scaleucd.o
+
clean:
/bin/rm -f *.o
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readnetgen.par
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readnetgen.par 2006-07-12 20:36:01 UTC (rev 4012)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readnetgen.par 2006-07-12 21:51:22 UTC (rev 4013)
@@ -0,0 +1,195 @@
+#
+# Parameter file for use with program readucd.
+#
+# **********************************************************************
+#
+# Coordinate parameters:
+#
+# cscale = coordinate scaling factor.
+# ix = coordinate number corresponding to x.
+# iy = coordinate number corresponding to y.
+# iz = coordinate number corresponding to z.
+#
+# The integer parameters may be used to transpose coordinates if a
+# different coordinate system is desired.
+#
+# cscale ix iy iz
+#-----------------------------------------------------------------------
+ 1.0d0 1 2 3
+#-----------------------------------------------------------------------
+#
+# Coordinate units:
+#
+# Give a string specifying the units used after scaling.
+#
+#-----------------------------------------------------------------------
+ 'km'
+#-----------------------------------------------------------------------
+#
+#
+# **********************************************************************
+#
+# Boundary condition information:
+# This information is used to determine the boundaries on which BC
+# will be applied, as well as the default BC.
+# Note: The boundaries will be determined after the mesh
+# coordinates have been scaled and permuted.
+#
+# nbc = number of boundary condition surfaces. Note that
+# each surface must be marked with a unique BC #.
+#
+# nbc
+#-----------------------------------------------------------------------
+ 5
+#-----------------------------------------------------------------------
+#
+# BC numbers and corresponding boundary conditions to be applied to
+# each boundary:
+# This input will consist of nbc lines, one for each mesh boundary.
+# The order is the same as for the mesh limits above.
+#
+# Each line will contain:
+# ibcode ibcx bcx ibcy bcy ibcz bcz acode
+#
+# where ibcode is the boundary condition code for a given surface,
+# each ibc value is the integer boundary condition code
+# (same as those for lithomop or tecton), and bc is the corresponding
+# boundary condition value.
+# The iauxcode value determines whether the actual BC value is read
+# from an auxiliary file or not. A value of zero indicates that
+# the specified bc value is the one to use. A value of one indicates
+# that the value is to come from the auxiliary file.
+#
+# ibcode ibcx(i) bcx(i) ibcy(i) bcy(i) ibcz(i) bcz(i) acode(i), i=1,nbc
+#-----------------------------------------------------------------------
+ 1 1 -999.0d0 1 -999.0d0 1 -999.0d0 1
+ 2 1 -999.0d0 1 -999.0d0 1 -999.0d0 1
+ 3 0 0.0d0 1 0.0d0 0 0.0d0 0
+ 4 1 -999.0d0 1 -999.0d0 1 -999.0d0 1
+ 5 1 -999.0d0 1 -999.0d0 1 -999.0d0 1
+#-----------------------------------------------------------------------
+#
+# Units used for boundary conditions:
+# You must specify the units to be used for displacement, velocity,
+# and force boundary conditions.
+#
+#-----------------------------------------------------------------------
+'m'
+# displacement_units
+'m/s'
+# velocity_units
+'newton'
+# force_units
+#-----------------------------------------------------------------------
+#
+# Connectivity order option:
+# This determines how the tetrahedral element connectivities are
+# ordered/determined.
+#
+# iconopt: 1 = assume connectivities are already in the proper
+# order (fastest, but frequently incorrect).
+# 2 = assume connectivities are in the wrong order. In
+# this case, the positions of nodes 2 and 3 are
+# reversed.
+# 3 = assume that the connectivities are in any order. In
+# this case, the element jacobian is first computed
+# using the original order. If a negative value is
+# obtained, the positions of nodes 2 and 3 are
+# reversed. This is probably the safest option.
+#
+# iconopt
+#-----------------------------------------------------------------------
+ 3
+#-----------------------------------------------------------------------
+#
+# Fault parameters:
+#
+# numflt = number of faults in the model
+# ibcflt = boundary condition code used to specify whether a node
+# lies on a fault
+#
+# numflt ibcflt
+#-----------------------------------------------------------------------
+ 1 7
+#-----------------------------------------------------------------------
+#
+# Fault specification:
+# Whether a node lies on a fault is determined by the attribute
+# information above. The fault to which a node belongs and the
+# side of the fault on which an element lies is determined by
+# material type. Thus, the two sides of a fault must always have
+# different material types, and if there is more than one fault,
+# a different material must lie along each side of a fault
+# boundary. For each fault, there are 3 lines of information:
+#
+# iftype(i), ifhist(i), nmatf(1,i), nmatf(2,i)
+# fsplit(j,1,i),j=1,nsd, ifmat(1,i,j),j=1,nmatf(1,i)
+# fsplit(j,2,i),j=1,nsd, ifmat(2,i,j),j=1,nmatf(2,i)
+#
+# The parameters are:
+#
+# iftype(i) = fault type (1=split node, 2=slippery node)
+# ifhist(i) = load history associated with the fault
+# (same type of specification as tecton/lithomop)
+# nmatf(1,i) = number of materials lying on side 1 of the fault
+# nmatf(2,i) = number of materials lying on side 2 of the fault
+# fsplit(j,1,i) = splitting parameter for each degree of freedom
+# on side 1 of the fault. If the fault is a split
+# node fault, this corresponds to an actual
+# dislocation value applied to this side of the
+# fault. If the fault is a slippery node fault,
+# the sign and magnitude are used to determine the
+# integer slip parameter (0 = no slip, -1 = slip
+# in negative local coordinate direction, 1 = slip
+# in positive local coordinate direction).
+# ifmat(1,i,j) = the list of material numbers lying on side 1 of
+# the fault.
+# fsplit(j,2,i) = same as fsplit(j,1,i) for side 2 of the fault.
+# ifmat(2,i,j) = same as ifmat(1,i,j) for side 2 of the fault.
+#
+# iftype(i), ifhist(i), nmatf(1,i), nmatf(2,i)
+# fsplit(j,1,i),j=1,nsd, ifmat(1,i,j),j=1,nmatf(1,i)
+# fsplit(j,2,i),j=1,nsd, ifmat(2,i,j),j=1,nmatf(2,i),i=1,numflt
+#-----------------------------------------------------------------------
+ 1 0 3 3
+ 0.3535d0 0.0d0 -0.3535d0 1 2 3
+ -0.3535d0 0.0d0 0.3535d0 4 5 6
+#-----------------------------------------------------------------------
+#
+# Winkler resoring force info:
+# In the present implementation, it is assumed that the Winkler BC
+# are being used to simulate gravity, although it would be easy to
+# adapt the BC to other uses. The user must first specify the
+# number of Winkler sets. If this number is zero, no additional
+# input is required.
+#
+# nwsets = number of Winkler restoring force sets.
+#
+# nwsets
+#-----------------------------------------------------------------------
+ 1
+#-----------------------------------------------------------------------
+#
+# Info describing each Winkler restoring force set. It is assumed
+# that Winkler BC for each set are applied in a single direction.
+# Enter one line for each Winkler BC set.
+#
+# The parameters are:
+#
+# iwcode(i) = The boundary condition code used to define this set.
+# iwdir(i) = The direction in which the force is to be applied.
+# 1 = X
+# 2 = Y
+# 3 = Z
+# modew(i) = Application mode code.
+# 0 = No restoring force (it would make no sense to
+# use this option).
+# 1 = Forces applied at all times.
+# -n = Forces applied according to load history n.
+# rhow(i) = Assumed density contrast.
+# gw(i) = Assumed graviational acceleration.
+#
+# iwcode(i) iwdir(i) modew(i) rhow(i) gw(i), i=1,nwsets
+#-----------------------------------------------------------------------
+ 6 3 1 3000.0 10.0
+#-----------------------------------------------------------------------
Property changes on: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readnetgen.par
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.f 2006-07-12 20:36:01 UTC (rev 4012)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.f 2006-07-12 21:51:22 UTC (rev 4013)
@@ -0,0 +1,552 @@
+ program readucd2
+c
+c... quick and dirty code to translate AVS UCD format to a format
+c suitable for lithomop or pylith. All boundaries for which BC
+c are to be applied should be flagged with a boundary condition
+c code. This also applies to any faults in the model.
+c
+c At present, code is only set up for linear tetrahedra, although
+c this would be easy to change.
+c
+ implicit none
+c
+c... parameters
+c
+ integer nsd,ndof,maxnodes,maxelmts,nen,maxflts,maxfnodes
+ integer maxbnds,maxfelems,maxnattr,maxeattr,ietyp,inf
+ parameter(nsd=3,ndof=3,maxnodes=1000000,maxelmts=5000000,nen=4,
+ & maxflts=6,maxfnodes=100000,maxfelems=100000,maxbnds=30,
+ & maxnattr=20,maxeattr=20,ietyp=5,inf=0)
+ integer kti,kto,kr,kw
+ parameter(kti=5,kto=6,kr=10,kw=11)
+ double precision eps
+ parameter(eps=1.0d-7)
+c
+c... local constants
+c
+ integer icflip(4)
+ data icflip/1,3,2,4/
+ integer kww,kwfb(maxflts),kwfc(maxflts)
+ data kww/12/
+ data kwfb/12,13,14,15,16,17/
+ data kwfc/18,19,20,21,22,23/
+ integer izero
+ data izero/0/
+ double precision zero
+ data zero/0.0d0/
+c
+c... parameters read from parameter file
+c
+ integer nbc,iconopt,numflt,ibfield
+ integer ibcode(maxbnds),ibc(nsd,maxbnds),iac(maxbnds),isn(nsd)
+ integer iftype(maxflts),iffield,ifcode(maxflts),ifefield(maxflts)
+ double precision bc(nsd,maxbnds),fsplit(nsd,2,maxflts)
+ double precision cscale
+ character cunits*20,dunits*20,vunits*20,funits*20
+c
+c... filenames
+c
+ character fileroot*200,pfile*200,ifile*200,nfile*200,cfile*200
+ character bcfile*200,afile*200
+ character fbcfile*200,fbccfile*200
+c
+c... parameters and variables read from UCD file
+c
+ integer numnp,numel,nnattr,neattr,nmattr
+ integer ien(nen,maxelmts),mat(maxelmts)
+ double precision x(nsd,maxnodes),attrn(maxnattr),attre(maxeattr)
+c
+c... values read from auxiliary file
+c
+ double precision bca(ndof,maxnodes)
+ double precision xa,ya,za
+c
+c... external routines
+c
+ double precision tetcmp
+ integer nnblnk,nchar
+ external nnblnk,nchar,tetcmp
+c
+c... local variables
+c
+ double precision det,sgn,xl(nsd,nen),xtmp(nsd),bct(3)
+ integer ifnodes(maxflts,maxfnodes),ifelems(2,maxflts,maxfelems)
+ integer nfnodes(maxflts),nfelems(2,maxflts)
+ integer ientmp(nen)
+ integer ifhist(maxflts),itmp(maxnattr),idir(nsd)
+ integer ibcnode(maxnodes)
+ integer i,j,k,l,n,jj,i1,i2,j1,j2,nenl,nsdl,ndofl,nfltnodes,nf
+ integer iattr,kk
+ integer elem,node,nflip
+ integer ibct
+ character cstring*14,etype*3,cfnum*2,descr*10
+ character dstring*21,vstring*17,fstring*14
+ data cstring/"coord_units = "/
+ data dstring/"displacement_units = "/
+ data vstring/"velocity_units = "/
+ data fstring/"force_units = "/
+ character stout*50
+ logical aux
+c
+ aux=.false.
+ nenl=nen
+ nsdl=nsd
+ ndofl=ndof
+ write(kto,*) "Enter root name for all files. Both input and"
+ write(kto,*) "output files will all have this prefix:"
+ read(kti,"(a200)") fileroot
+ i1=nnblnk(fileroot)
+ i2=nchar(fileroot)
+ call ifill(ibcnode,izero,maxnodes)
+ call ifill(nfnodes,izero,maxflts)
+ call ifill(nfelems,izero,2*maxflts)
+c
+c... read parameter file
+c
+ pfile=fileroot(i1:i2)//".par"
+ open(file=pfile,unit=kr,status="old")
+c
+c... coordinate scaling factor and units
+c
+ call pskip(kr)
+ read(kr,*) cscale,(idir(i),i=1,nsd)
+ if(cscale.eq.0.0d0) cscale=1.0d0
+ call pskip(kr)
+ read(kr,*) cunits
+c
+c... bc info
+c
+ call pskip(kr)
+ read(kr,*) nbc,ibfield
+ do i=1,nbc
+ call pskip(kr)
+ read(kr,*) ibcode(i),(ibc(j,ibcode(i)),bc(j,ibcode(i)),j=1,nsd),
+ & iac(ibcode(i))
+ if(iac(ibcode(i)).ne.0) aux=.true.
+ end do
+ call pskip(kr)
+ read(kr,*) dunits
+ call pskip(kr)
+ read(kr,*) vunits
+ call pskip(kr)
+ read(kr,*) funits
+c
+c... connectivity order option
+c
+ call pskip(kr)
+ read(kr,*) iconopt
+c
+c... fault definitions
+c
+ call pskip(kr)
+ read(kr,*) numflt,iffield
+ do i=1,numflt
+ call pskip(kr)
+ read(kr,*) iftype(i),ifhist(i),ifcode(i),ifefield(i)
+ call pskip(kr)
+ read(kr,*) (fsplit(j,1,i),j=1,nsd)
+ call pskip(kr)
+ read(kr,*) (fsplit(j,2,i),j=1,nsd)
+ end do
+ close(kr)
+c
+c... read and output nodal coordinates
+c
+ ifile=fileroot(i1:i2)//".inp"
+ open(file=ifile,unit=kr,status="old")
+ nfile=fileroot(i1:i2)//".coord"
+ open(file=nfile,unit=kw,status="new")
+ stout=cstring//cunits
+ write(kw,"(a50)") stout
+ call pskip(kr)
+ read(kr,*) numnp,numel,nnattr,neattr,nmattr
+ do i=1,numnp
+ read(kr,*) n,(xtmp(j),j=1,nsd)
+ do j=1,nsd
+ x(j,i)=cscale*xtmp(idir(j))
+ end do
+ write(kw,"(i7,3(2x,1pe15.8))") i,(x(j,i),j=1,nsd)
+ end do
+ close(kw)
+c
+c... read and output connectivity info
+c
+ cfile=fileroot(i1:i2)//".connect"
+ open(file=cfile,unit=kw,status="new")
+ nflip=0
+ if(iconopt.eq.1) then
+ do i=1,numel
+ read(kr,*) n,mat(i),etype,(ien(j,i),j=1,nen)
+ write(kw,"(i7,3i4,4i7)") i,ietyp,mat(i),inf,(ien(j,i),j=1,nen)
+ end do
+ else if(iconopt.eq.2) then
+ do i=1,numel
+ read(kr,*) n,mat(i),etype,(ientmp(j),j=1,nen)
+ do j=1,nen
+ ien(j,i)=ientmp(icflip(j))
+ end do
+ write(kw,"(i7,3i4,4i7)") i,ietyp,mat(i),inf,(ien(j,i),j=1,nen)
+ end do
+ else if(iconopt.eq.3) then
+ do i=1,numel
+ read(kr,*) n,mat(i),etype,(ientmp(j),j=1,nen)
+ call lcoord(ientmp,x,xl,nsdl,nenl,numnp)
+ det=tetcmp(xl)
+ if(det.lt.0.0d0) then
+ nflip=nflip+1
+ do j=1,nen
+ ien(j,i)=ientmp(icflip(j))
+ end do
+ else
+ do j=1,nen
+ ien(j,i)=ientmp(j)
+ end do
+ end if
+ write(kw,"(i7,3i4,4i7)") i,ietyp,mat(i),inf,(ien(j,i),j=1,nen)
+ end do
+ end if
+ close(kw)
+c
+c... read nodal attributes to determine which nodes are associated
+c with each fault and boundary condition code.
+c
+ read(kr,*) nf,(itmp(i),i=1,nf)
+ do i=1,nf
+ read(kr,*) descr
+ end do
+ nfltnodes=0
+ do i=1,numnp
+ read(kr,*) n,(attrn(j),j=1,nnattr)
+ iattr=nint(attrn(ibfield))
+ do j=1,nbc
+ if(iattr.eq.ibcode(j)) ibcnode(i)=iattr
+ end do
+c
+ iattr=nint(attrn(iffield))
+ do j=1,numflt
+ if(iattr.eq.ifcode(j)) then
+ nfnodes(j)=nfnodes(j)+1
+ nfltnodes=nfltnodes+1
+ ifnodes(j,nfnodes(j))=i
+ end if
+ end do
+ end do
+c
+c... read element attributes to determine which elements are adjacent
+c to each fault
+c
+ read(kr,*) nf,(itmp(i),i=1,nf)
+ do i=1,nf
+ read(kr,*) descr
+ end do
+ do i=1,numel
+ read(kr,*) n,(attre(j),j=1,neattr)
+ do j=1,numflt
+ iattr=nint(attre(ifefield(j)))
+ if(iattr.eq.1) then
+ nfelems(1,j)=nfelems(1,j)+1
+ ifelems(1,j,nfelems(1,j))=n
+ else if(iattr.eq.-1) then
+ nfelems(2,j)=nfelems(2,j)+1
+ ifelems(2,j,nfelems(2,j))=n
+ end if
+ end do
+ end do
+ close(kr)
+c
+c... if auxiliary file is being used, read BC from it
+c
+ if(aux) then
+ afile=fileroot(i1:i2)//".aux"
+ open(file=afile,unit=kr,status="old")
+ do i=1,numnp
+ read(kr,*) xa,ya,za,(bca(j,i),j=1,ndofl)
+ end do
+ close(kr)
+ end if
+c
+c... output BC and BC coordinates
+c
+ bcfile=fileroot(i1:i2)//".bc"
+ open(file=bcfile,unit=kw,status="new")
+ stout=dstring//dunits
+ write(kw,"(a50)") stout
+ stout=vstring//vunits
+ write(kw,"(a50)") stout
+ stout=fstring//funits
+ write(kw,"(a50)") stout
+ do i=1,numnp
+ ibct=ibcnode(i)
+ if(ibct.ne.0) then
+ if(iac(ibct).eq.0) then
+ do j=1,3
+ bct(j)=bc(j,ibct)
+ end do
+ else
+ do j=1,3
+ bct(j)=bca(j,i)
+ end do
+ end if
+ write(kw,"(i7,3(2x,i3),3(2x,1pe15.8))") i,
+ & (ibc(k,ibct),k=1,nsd),(bct(k),k=1,nsd)
+ end if
+ end do
+ close(kw)
+c
+c... output fault info after determining which fault each node lies on
+c
+ do i=1,numflt
+ write(cfnum,"(i2)") i
+ j1=nnblnk(cfnum)
+ j2=nchar(cfnum)
+ fbcfile=fileroot(i1:i2)//"."//cfnum(j1:j2)//".fbc"
+ fbccfile=fileroot(i1:i2)//"."//cfnum(j1:j2)//".fcoord"
+ open(file=fbcfile,unit=kwfb(i),status="new")
+ open(file=fbccfile,unit=kwfc(i),status="new")
+ end do
+c
+c... Loop over positive sides of fault, then negative.
+c
+ do kk=1,2
+ do i=1,numflt
+ do j=1,nsd
+ sgn=sign(1.0d0,fsplit(j,kk,i))
+ isn(j)=nint(sgn)
+ if(fsplit(j,kk,i).eq.0.0d0) isn(j)=0
+ end do
+c
+c... First loop over faulted elements, then nodes attached to each element.
+c
+ do j=1,nfelems(1,i)
+ elem=ifelems(1,i,j)
+ do l=1,nen
+ do k=1,nfnodes(i)
+ node=ifnodes(i,k)
+ if(ien(l,elem).eq.node) then
+ if(iftype(i).eq.1) then
+ write(kwfb(i),"(2i7,i4,3(2x,1pe15.8))") elem,
+ & node,ifhist(i),(fsplit(jj,kk,i),jj=1,nsd)
+ else
+ write(kwfb(i),"(2i7,3i4)") elem,node,
+ & (isn(jj),jj=1,nsd)
+ end if
+ write(kwfc(i),"(3i7,3(2x,1pe15.8))")
+ & elem,node,kk,(x(jj,node),jj=1,nsd)
+ end if
+ end do
+ end do
+ end do
+ end do
+ end do
+ do i=1,numflt
+ close(kwfb(i))
+ close(kwfc(i))
+ end do
+ write(kto,700) nflip
+700 format("Number of connectivities flipped: ",i7)
+ stop
+ end
+c
+c
+ subroutine ifill(iarr,ival,nlen)
+c
+c... subroutine to fill an integer array with a given value
+c
+ implicit none
+ integer ival,nlen
+ integer iarr(nlen)
+ integer i
+ do i=1,nlen
+ iarr(i)=ival
+ end do
+ return
+ end
+c
+c
+ subroutine lcoord(ien,x,xl,nsd,nen,numnp)
+c
+c... subroutine to localize element coordinates
+c
+ implicit none
+ integer nsd,nen,numnp
+ integer ien(nen)
+ double precision x(nsd,numnp),xl(nsd,nen)
+c
+ integer i,j,ii
+c
+ do i=1,nen
+ ii=ien(i)
+ do j=1,nsd
+ xl(j,i)=x(j,ii)
+ end do
+ end do
+ return
+ end
+c
+c
+ function nchar(string)
+c
+c... determines the minimum nonblank length of a string
+c
+ implicit none
+c
+c... parameter definitions
+c
+ character blank*1
+ parameter(blank=' ')
+c
+c... function arguments
+c
+ integer nchar
+ character*(*) string
+c
+c... intrinsic functions
+c
+ intrinsic len
+c
+c... local variables
+c
+ integer nmax,i,itest
+c
+ nmax=len(string)
+ nchar=0
+ do i=1,nmax
+ itest=nmax-i+1
+ if(string(itest:itest).ne.blank) then
+ nchar=itest
+ return
+ end if
+ end do
+ return
+ end
+c
+c
+ function nnblnk(string)
+c
+c determines the position of the first nonblank entry
+c of a string (returns 1 if the first character is
+c not blank)
+c
+ implicit none
+c
+c... parameter definitions
+c
+ character blank*1
+ parameter(blank=' ')
+c
+c... function arguments
+c
+ integer nnblnk
+ character*(*) string
+c
+c... intrinsic functions
+c
+ intrinsic len
+c
+c... local variables
+c
+ integer nmax,i
+c
+ nmax=len(string)
+ nnblnk=nmax
+ do i=1,nmax
+ if(string(i:i).ne.blank) then
+ nnblnk=i
+ return
+ end if
+ end do
+ return
+ end
+c
+c
+ subroutine pskip(iunit)
+c
+c routine to skip lines beginning with the string # and blank
+c lines.
+c this routine ignores leading blanks before the key string.
+c
+c
+ implicit none
+c
+c... subroutine arguments
+c
+ integer iunit
+c
+c... local constants
+c
+ character leader*1
+ data leader/'#'/
+c
+c... intrinsic functions
+c
+ intrinsic index
+c
+c... external functions
+c
+ integer nchar,nnblnk
+ external nchar,nnblnk
+c
+c... local variables
+c
+ integer inblnk
+ character string*80
+c
+ 10 continue
+ read(iunit,"(a80)",end=20) string
+ if(nchar(string).eq.0) goto 10
+ inblnk=nnblnk(string)
+ if(index(string,leader).eq.inblnk) goto 10
+ backspace(iunit)
+ 20 continue
+ return
+ end
+c
+c
+ function tetcmp(xl)
+c
+c... program to compute determinant and equivalent tetrahedral volume
+c for a given 4x4 matrix.
+c
+ implicit none
+ double precision tetcmp,xl(3,4)
+ double precision a(4,4),cof(3,3)
+ data a(1,1),a(1,2),a(1,3),a(1,4)/1.0d0,1.0d0,1.0d0,1.0d0/
+ integer i,j,ii
+ double precision det3,onem,sgn
+ external det3
+ integer ind(3,4)
+ data ind/2,3,4,1,3,4,1,2,4,1,2,3/
+ do i=2,4
+ do j=1,4
+ a(i,j)=xl(i-1,j)
+ end do
+ end do
+ onem=-1.0d0
+ sgn=onem
+ tetcmp=0.0d0
+ do i=1,4
+ sgn=sgn*onem
+ do ii=1,3
+ do j=2,4
+ cof(j-1,ii)=a(j,ind(ii,i))
+ end do
+ end do
+ tetcmp=tetcmp+sgn*a(1,i)*det3(cof)
+cdebug write(6,*) "i,a,tetcmp:",i,(a(j,i),j=1,4)
+ end do
+ return
+ end
+c
+c
+ function det3(x)
+c
+c... function to compute determinant of a 3x3 matrix
+c
+ implicit none
+ double precision det3
+ double precision x(3,3)
+ det3=x(1,1)*x(2,2)*x(3,3)-x(1,1)*x(2,3)*x(3,2)+
+ & x(1,2)*x(2,3)*x(3,1)-x(1,2)*x(2,1)*x(3,3)+
+ & x(1,3)*x(2,1)*x(3,2)-x(1,3)*x(2,2)*x(3,1)
+ return
+ end
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.par
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.par 2006-07-12 20:36:01 UTC (rev 4012)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.par 2006-07-12 21:51:22 UTC (rev 4013)
@@ -0,0 +1,164 @@
+#
+# Parameter file for use with program readucd.
+#
+# **********************************************************************
+#
+# Coordinate parameters:
+#
+# cscale = coordinate scaling factor.
+# ix = coordinate number corresponding to x.
+# iy = coordinate number corresponding to y.
+# iz = coordinate number corresponding to z.
+#
+# The integer parameters may be used to transpose coordinates if a
+# different coordinate system is desired.
+#
+# cscale ix iy iz
+#-----------------------------------------------------------------------
+ 1.0d0 1 2 3
+#-----------------------------------------------------------------------
+#
+# Coordinate units:
+#
+# Give a string specifying the units used after scaling.
+#
+#-----------------------------------------------------------------------
+ 'km'
+#-----------------------------------------------------------------------
+#
+#
+# **********************************************************************
+#
+# Boundary condition information:
+# This information is used to determine the boundaries on which BC
+# will be applied, as well as the default BC.
+#
+# nbc = number of boundary condition surfaces. Note that
+# each surface must be marked with a unique BC #.
+# ibfield= field number from which BC codes are read. This must be
+# a node-based field (attribute) in the UCD file.
+#
+# nbc ibfield
+#-----------------------------------------------------------------------
+ 5 5
+#-----------------------------------------------------------------------
+#
+# BC numbers and corresponding boundary conditions to be applied to
+# each boundary:
+# This input will consist of nbc lines, one for each mesh boundary.
+# The order is the same as for the mesh limits above.
+#
+# Each line will contain:
+# ibcode ibcx bcx ibcy bcy ibcz bcz acode
+#
+# where ibcode is the boundary condition code for a given surface,
+# each ibc value is the integer boundary condition code
+# (same as those for lithomop or tecton), and bc is the corresponding
+# boundary condition value.
+# The iauxcode value determines whether the actual BC value is read
+# from an auxiliary file or not. A value of zero indicates that
+# the specified bc value is the one to use. A value of one indicates
+# that the value is to come from the auxiliary file.
+#
+# ibcode ibcx(i) bcx(i) ibcy(i) bcy(i) ibcz(i) bcz(i) acode(i), i=1,nbc
+#-----------------------------------------------------------------------
+ 1 0 0.0d0 1 0.0d0 0 0.0d0 0
+ 2 1 -999.0d0 1 -999.0d0 1 -999.0d0 1
+ 3 1 -999.0d0 1 -999.0d0 1 -999.0d0 1
+ 4 1 -999.0d0 1 -999.0d0 1 -999.0d0 1
+ 5 1 -999.0d0 1 -999.0d0 1 -999.0d0 1
+# The numbering above corresponds to front, back, left, right, and
+# bottom in the UCD file.
+#-----------------------------------------------------------------------
+#
+# Units used for boundary conditions:
+# You must specify the units to be used for displacement, velocity,
+# and force boundary conditions.
+#
+#-----------------------------------------------------------------------
+'m'
+# displacement_units
+'m/s'
+# velocity_units
+'newton'
+# force_units
+#-----------------------------------------------------------------------
+#
+# Connectivity order option:
+# This determines how the tetrahedral element connectivities are
+# ordered/determined.
+#
+# iconopt: 1 = assume connectivities are already in the proper
+# order (fastest, but frequently incorrect).
+# 2 = assume connectivities are in the wrong order. In
+# this case, the positions of nodes 2 and 3 are
+# reversed.
+# 3 = assume that the connectivities are in any order. In
+# this case, the element jacobian is first computed
+# using the original order. If a negative value is
+# obtained, the positions of nodes 2 and 3 are
+# reversed. This is probably the safest option.
+#
+# iconopt
+#-----------------------------------------------------------------------
+ 3
+#-----------------------------------------------------------------------
+#
+# Fault parameters:
+#
+# numflt = number of faults in the model
+# iffield= field number from which fault info is determined. This
+# must be a node-based field (attribute) in the UCD file.
+#
+# numflt iffield
+#-----------------------------------------------------------------------
+ 1 6
+#-----------------------------------------------------------------------
+#
+# Fault specification:
+# Whether a node lies on a fault is determined by the attribute
+# information above. The fault to which a node belongs (if any) is
+# determined by the value of the iffield entry for each node.
+# The fault to which an element belongs and the side (+ or -) is
+# determined by the entry in ifefield as described below. Each
+# fault has its element info determined by a separate element-based
+# field. A zero entry in this field indicates the element is not
+# attached to a fault, a -1 indicates the element is on the negative
+# side of the fault, and a +1 indicates the element is on the
+# positive side of the fault.
+# For each fault, there are 3 lines of information:
+#
+# iftype(i), ifhist(i), ifcode(i), ifefield(i)
+# fsplit(j,1,i),j=1,nsd
+# fsplit(j,2,i),j=1,nsd
+#
+# The parameters are:
+#
+# iftype(i) = fault type (1=split node, 2=slippery node)
+# ifhist(i) = load history associated with the fault
+# (same type of specification as tecton/lithomop)
+# ifcode(i) = the nodal code value corresponding to this fault,
+# read from iffield above. Nodes with iffield
+# values corresponding to this value are on the
+# fault.
+# ifefield(i) = the element field value corresponding to this
+# fault.
+# fsplit(j,1,i) = splitting parameter for each degree of freedom
+# on side 1 of the fault. If the fault is a split
+# node fault, this corresponds to an actual
+# dislocation value applied to this side of the
+# fault. If the fault is a slippery node fault,
+# the sign and magnitude are used to determine the
+# integer slip parameter (0 = no slip, -1 = slip
+# in negative local coordinate direction, 1 = slip
+# in positive local coordinate direction).
+# fsplit(j,2,i) = same as fsplit(j,1,i) for side 2 of the fault.
+#
+# iftype(i), ifhist(i), ifcode(i), ifefield(i)
+# fsplit(j,1,i),j=1,nsd
+# fsplit(j,2,i),j=1,nsd
+#-----------------------------------------------------------------------
+ 1 0 1 1
+ 0.3535d0 0.0d0 -0.3535d0
+ -0.3535d0 0.0d0 0.3535d0
+#-----------------------------------------------------------------------
Property changes on: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/readucd2.par
___________________________________________________________________
Name: svn:executable
+ *
Added: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/scaleucd.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/scaleucd.f 2006-07-12 20:36:01 UTC (rev 4012)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/scaleucd.f 2006-07-12 21:51:22 UTC (rev 4013)
@@ -0,0 +1,141 @@
+ program scaleucd
+c
+c... Program to scale the coordinates produced by netgen and output
+c them to a separate file.
+c
+ implicit none
+c
+c... filenames and unit numbers
+c
+ integer kti,kto,kr,kw
+ common/units/kti,kto,kr,kw
+ character ifile*200,ofile*200
+c
+c... local variables
+c
+ integer numnp,numel,nnattr,neattr,nmattr,n,i,j
+ double precision x(3),scale
+c
+c... get filenames and scale factor
+c
+ call files(scale,ifile,ofile)
+ open(kr,file=ifile,status="old")
+ open(kw,file=ofile,status="new")
+c
+c... get number of nodes and loop over them
+c
+ read(kr,*) numnp,numel,nnattr,neattr,nmattr
+ do i=1,numnp
+ read(kr,*) n,(x(j),j=1,3)
+ write(kw,"(3(2x,1pe15.8))") (scale*x(j),j=1,3)
+ end do
+ close(kr)
+ close(kw)
+ stop
+ end
+c
+c
+ subroutine files(scale,ifile,ofile)
+c
+c... subroutine to get filenames and scale factor from command-line
+c
+ implicit none
+ double precision scale
+ character ifile*(*),ofile*(*)
+c
+c... intrinsic functions
+c
+ intrinsic iargc,index
+c
+c... external functions
+c
+ integer nchar
+ external nchar
+c
+c... unit numbers
+c
+ integer kti,kto,kr,kw
+ common/units/kti,kto,kr,kw
+c
+c... local variables
+c
+ integer nargs,i,j
+ character string*2000
+ logical fflag(3)
+c
+ kti=5
+ kto=6
+ kr=10
+ kw=11
+c
+ do i=1,3
+ fflag(i)=.false.
+ end do
+c
+ nargs=iargc()
+ do i=1,nargs
+ call getarg(i,string)
+ if(index(string,'i=').ne.0) then
+ j=nchar(string)
+ ifile=string(3:j)
+ fflag(1)=.true.
+ else if(index(string,'o=').ne.0) then
+ j=nchar(string)
+ ofile=string(3:j)
+ fflag(2)=.true.
+ else if(index(string,'s=').ne.0) then
+ j=nchar(string)
+ read(string(3:j),*) scale
+ fflag(3)=.true.
+ end if
+ end do
+c
+ do i=1,3
+ if(.not.fflag(i)) then
+ write(kto,800)
+ stop
+ end if
+ end do
+c
+ 800 format("Usage:",/,
+ & " scaleucd i=<input_file> o=<output_file>",/,
+ & " s=<scale_factor>")
+ return
+ end
+c
+c
+ function nchar(string)
+c
+c... determines the minimum nonblank length of a string
+c
+ implicit none
+c
+c... parameter definitions
+c
+ character blank*1
+ parameter(blank=' ')
+c
+c... function arguments
+c
+ integer nchar
+ character*(*) string
+c
+c... intrinsic functions
+c
+ intrinsic len
+c
+c... local variables
+c
+ integer nmax,i,itest
+c
+ nmax=len(string)
+ nchar=0
+ do i=1,nmax
+ itest=nmax-i+1
+ if(string(itest:itest).ne.blank) then
+ nchar=itest
+ return
+ end if
+ end do
+ return
+ end
Property changes on: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/scaleucd.f
___________________________________________________________________
Name: svn:executable
+ *
More information about the cig-commits
mailing list