[cig-commits] r4227 -
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils
willic3 at geodynamics.org
willic3 at geodynamics.org
Fri Aug 4 13:51:39 PDT 2006
Author: willic3
Date: 2006-08-04 13:51:39 -0700 (Fri, 04 Aug 2006)
New Revision: 4227
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/makeucd.f
Log:
Initial version of code to turn UCD pieces from LithoMop/PyLith into
complete UCD files. All behavior is controlled by command-line
arguments including options for interpolating state variables to nodes
and averaging Gauss point values to element centroids.
The code definitely needs more testing. So far, it produces a valid
UCD file and the results look reasonable, but no quantitative testing
has been done.
Likely future enhancements include the ability to deal with binary UCD
files.
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/makeucd.f
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/makeucd.f 2006-08-04 19:18:23 UTC (rev 4226)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/utils/makeucd.f 2006-08-04 20:51:39 UTC (rev 4227)
@@ -25,8 +25,8 @@
c... input/output and logical flag parameters
c
integer kti,kto,km,kp,kg,ko
- common/kti,kto,km,kp,kg,ko/
- logical passigned,gassigned,ninterp,gavg,nout,cout
+ common/units/kti,kto,km,kp,kg,ko
+ logical passigned,gassigned,ninterp,nout,cout
c
c... parameters and variables read/deduced from mesh UCD file
c
@@ -37,7 +37,7 @@
c... parameters and variables read/deduced from nodal values file
c
integer nnattr,isnattr(maxattr),ntnattr
- double precision vnodes(maxattr*naxnodes)
+ double precision vnodes(maxattr*maxnodes)
character nvnames(maxattr)*30
c
c... parameters and variables read/deduced from element values file
@@ -50,7 +50,8 @@
c
integer nnattrout,neattrout,ntnattrout,nteattrout
integer isnattrout(maxoattr),iseattrout(maxoattr)
- double precision vnout(maxoattr*maxnodes),veout(maxoattr*maxelmts)
+ double precision vnout(maxoattr*maxnodes),vavg(maxoattr)
+ character nvnamesout(maxoattr)*30,evnamesout(maxoattr)*30
c
c... external routines
c
@@ -68,16 +69,20 @@
c
c... local variables
c
- integer idum,i,nen,ngauss,ind
+ integer nsdl
+ integer idum,i,j,ind,n,ietype
integer ientmp(maxnen)
character tstring*500,etype*3
+ double precision smop(maxnodes)
double precision sh((nsd+1)*maxnen*maxgauss)
- double precision gauss((nsd+1)*maxgauss)
+ double precision gauss((nsd+1)*maxgauss),xl(nsd*maxnen)
c
+ nsdl=nsd
+c
c... get command-line arguments, open files, and assign logical
c variables
c
- call files(passigned,gassigned,ninterp,gavg,nout,cout)
+ call files(passigned,gassigned,ninterp,nout,cout)
c
c... read mesh UCD file
c
@@ -85,11 +90,14 @@
do i=1,numnp
read(km,*) n,(x(j,n),j=1,nsd)
end do
- read(km,*) tstring
+ read(km,"(a500)") tstring
+ ietype=0
if(index(tstring,"hex").ne.0) then
+ ietype=1
nen=8
ngauss=8
else if(index(tstring,"tet").ne.0) then
+ ietype=2
nen=4
ngauss=1
else
@@ -102,7 +110,7 @@
do i=1,numel
read(km,*) n,mat(n),etype,(ientmp(j),j=1,nen)
do j=1,nen
- ien(j+ind)=ientmp(indp(j))
+ ien(j+ind)=ientmp(indp(j,ietype))
end do
ind=ind+nen
end do
@@ -153,9 +161,28 @@
nteattrout=0
if(nout) then
nnattrout=nnattr
- if(ninterp) nnattrout=nnattrout+neattr
+ ntnattrout=ntnattr
+ do i=1,nnattr
+ isnattrout(i)=isnattr(i)
+ nvnamesout(i)=nvnames(i)
+ end do
+ if(ninterp) then
+ do i=1,neattr
+ isnattrout(i+nnattrout)=iseattr(i)
+ nvnamesout(i+nnattrout)=evnames(i)
+ end do
+ nnattrout=nnattrout+neattr
+ ntnattrout=ntnattrout+nteattr
+ end if
end if
- if(gavg) neattrout=neattr
+ if(cout) then
+ neattrout=neattr
+ nteattrout=nteattr
+ do i=1,neattr
+ iseattrout(i)=iseattr(i)
+ evnamesout(i)=evnames(i)
+ end do
+ end if
idum=0
write(ko,"(5i9)") numnp,numel,nnattrout,neattrout,idum
c
@@ -170,228 +197,91 @@
ientmp(j)=ien(ind+j)
end do
write(ko,"(2i8,2x,a3,10i9)") i,mat(i),etype,
- & (ientmp(indu(j)),j=1,nen)
+ & (ientmp(indu(j,ietype)),j=1,nen)
ind=ind+nen
end do
-
- nenl=nen
- nsdl=nsd
- 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)
c
-c... read parameter file
+c... output nodal values if requested, interpolating Gauss point
+c values to nodes if desired.
c
- pfile=fileroot(i1:i2)//".par"
- open(file=pfile,unit=kr,status="old")
+ if(nout) call nodout(
+ & x,ien,numnp,numel,nsdl,nen,ngauss,ninterp,
+ & vnodes,isnattr,nnattr,ntnattr,
+ & velemt,iseattr,neattr,nteattr,
+ & vnout,isnattrout,nnattrout,ntnattrout,nvnamesout,
+ & smop,sh,gauss,xl)
c
-c... coordinate scaling factor
+c... output values at element centroids if requested
c
- call pskip(kr)
- read(kr,*) cscale,(idir(i),i=1,nsd)
- if(cscale.eq.0.0d0) cscale=1.0d0
+ if(cout) call centout(
+ & numel,ngauss,
+ & velemt,neattr,nteattr,
+ & vavg,iseattrout,neattrout,nteattrout,evnamesout)
c
-c... box dimensions and bc
+ close(ko)
c
- call pskip(kr)
- read(kr,*) ((xlim(j,i),j=1,2),i=1,nsd)
- do i=1,nsides
- call pskip(kr)
- read(kr,*) (ibc(j,i),bc(j,i),j=1,nsd)
- end do
-c... connectivity order option
+ stop
+ end
c
- call pskip(kr)
- read(kr,*) iconopt
c
-c... fault definitions
+ subroutine centout(
+ & numel,ngauss,
+ & velemt,neattr,nteattr,
+ & vavg,iseattrout,neattrout,nteattrout,evnamesout)
c
- call pskip(kr)
- read(kr,*) numflt,ifattrn,ifattrv
- do i=1,numflt
- call pskip(kr)
- read(kr,*) iftype(i),ifhist(i),nmatf(1,i),nmatf(2,i)
- call pskip(kr)
- read(kr,*) (fsplit(j,1,i),j=1,nsd),(ifmat(1,i,j),j=1,nmatf(1,i))
- call pskip(kr)
- read(kr,*) (fsplit(j,2,i),j=1,nsd),(ifmat(2,i,j),j=1,nmatf(2,i))
- end do
- close(kr)
+c... routine to output values at element centroids
c
-c... read and output nodal coordinates and bc info
+ implicit none
+ integer numel,ngauss
+ integer neattr,nteattr,neattrout,nteattrout
+ integer iseattrout(neattrout)
+ double precision velemt(nteattr,ngauss,numel)
+ double precision vavg(nteattrout)
+ character evnamesout(neattrout)*(*)
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")
- do i=1,nsides
- j1=nnblnk(cside(i))
- j2=nchar(cside(i))
- bcfile=fileroot(i1:i2)//"."//cside(i)(j1:j2)//".bc"
- bccfile=fileroot(i1:i2)//"."//cside(i)(j1:j2)//".coord"
- open(file=bcfile,unit=kwb(i),status="new")
- open(file=bccfile,unit=kwc(i),status="new")
- end do
- write(kw,"(a15)") cstring
- 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)
- jj=0
- do j=1,nsd
- do k=1,2
- jj=jj+1
- dr=abs(x(j,i)-xlim(k,j))
- if(dr.lt.eps) then
- write(kwb(jj),"(i7,3i4,3(2x,1pe15.8))") i,
- & (ibc(l,jj),l=1,nsd),(bc(l,jj),l=1,nsd)
- write(kwc(jj),"(i7,3(2x,1pe15.8))") i,
- & (x(l,i),l=1,nsd)
- end if
- end do
- end do
- end do
- close(kw)
- do i=1,nsides
- close(kwb(i))
- close(kwc(i))
- end do
+c... i/o units
c
-c... read and output connectivity info
+ integer kti,kto,km,kp,kg,ko
+ common/units/kti,kto,km,kp,kg,ko
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)
-cdebug write(kto,*) "i,det:",i,det
- 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... local variables
c
-c... read nodal attributes to determine which nodes are associated
-c with each fault
+ integer i,j,k
+ double precision wt
c
- read(kr,*) nf,(itmp(i),i=1,nf)
- do i=1,nf
- read(kr,*) descr
+c... output cell values header
+c
+ write(ko,"(30i5)") neattrout,(iseattrout(i),i=1,neattrout)
+ do i=1,neattrout
+ write(ko,"(a30)") evnamesout(i)
end do
- nfltnodes=0
- do i=1,numnp
- read(kr,*) n,(fattr(j),j=1,nnattr)
- iattr=nint(fattr(ifattrn))
- if(iattr.eq.ifattrv) then
- nfltnodes=nfltnodes+1
- inodef(nfltnodes)=i
- nelsf(nfltnodes)=0
- end if
- end do
c
-c... determine which elements contain each node on the faults
+c... compute values averaged over an element and output them
c
+ wt=1.0d0/dble(ngauss)
do i=1,numel
- do j=1,nfltnodes
- do k=1,nen
- if(ien(k,i).eq.inodef(j)) then
- nelsf(j)=nelsf(j)+1
- iadjf(j,nelsf(j))=i
- end if
+ call fill(vavg,0.0d0,nteattrout)
+ do j=1,ngauss
+ do k=1,nteattr
+ vavg(k)=vavg(k)+wt*velemt(k,j,i)
end do
end do
+ write(ko,"(i7,30(2x,1pe15.8))") i,(vavg(j),j=1,nteattrout)
end do
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... First find all elements on one side of the fault, then the other
-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)
- end do
- do j=1,nfltnodes
- do k=1,nelsf(j)
- do l=1,nmatf(kk,i)
- iel=iadjf(j,k)
- if(mat(iel).eq.ifmat(kk,i,l)) then
- if(iftype(i).eq.1) then
- write(kwfb(i),"(2i7,i4,3(2x,1pe15.8))") iel,
- & inodef(j),ifhist(i),(fsplit(jj,kk,i),jj=1,nsd)
- else
- write(kwfb(i),"(2i7,3i4)") iel,inodef(j),
- & (isn(jj),jj=1,nsd)
- end if
- write(kwfc(i),"(3i7,3(2x,1pe15.8))")
- & iel,inodef(j),kk,(x(jj,inodef(j)),jj=1,nsd)
- end if
- end do
- end do
- end do
- end do
- end do
- do i=1,nsides
- close(kwfb(i))
- close(kwfc(i))
- end do
- write(kto,700) nflip
-700 format("Number of connectivities flipped: ",i7)
- stop
+ return
end
c
c
- subroutine files(passigned,gassigned,ninterp,gavg,nout,cout)
+ subroutine files(passigned,gassigned,ninterp,nout,cout)
c
c... subroutine to set up i/o and run options
c
implicit none
- logical ninterp,gavg,nout,cout
+ logical ninterp,nout,cout
c
integer kti,kto,km,kp,kg,ko
- common/kti,kto,km,kp,kg,ko/
+ common/units/kti,kto,km,kp,kg,ko
c
intrinsic index,iargc
c
@@ -414,7 +304,6 @@
c... logical flags
c
ninterp=.true.
- gavg=.false.
nout=.true.
cout=.false.
massigned=.false.
@@ -452,21 +341,10 @@
ninterp=.true.
else
write(kto,800)
- end
+ stop
end if
- else if(index(string,'a=').ne.0) then
+ else if(index(string,'nf=').ne.0) then
j=lnblnk(string)
- read(string(3:j),*) ival
- if(ival.eq.0) then
- gavg=.false.
- else if(ival.eq.1) then
- gavg=.true.
- else
- write(kto,800)
- end
- end if
- else if(index(string,'no=').ne.0) then
- j=lnblnk(string)
read(string(4:j),*) ival
if(ival.eq.0) then
nout=.false.
@@ -474,9 +352,9 @@
nout=.true.
else
write(kto,800)
- end
+ stop
end if
- else if(index(string,'co=').ne.0) then
+ else if(index(string,'cf=').ne.0) then
j=lnblnk(string)
read(string(4:j),*) ival
if(ival.eq.0) then
@@ -485,7 +363,7 @@
cout=.true.
else
write(kto,800)
- end
+ stop
end if
end if
end do
@@ -493,7 +371,6 @@
c... adjust flags for inconsistencies
c
if(.not.nout) ninterp=.false.
- if(.not.cout) gavg=.false.
c
c... open files
c
@@ -515,19 +392,77 @@
800 format("Usage:",/,
& "makeucd m=<mesh_file> o=<output_file> [p=<nodal_value_file>]",/,
& "[g=<gauss_value_file>] [n=<gauss_interpolation_option]",/,
- & "[a=<gauss_averaging_option>] [no=<nodal_output_option>]",/,
- & "[co=<element_output_option>]",//,
+ & "[nf=<nodal_output_flag>] [cf=<element_output_flag>]",//,
& "All options have values of 0 (no) or 1 (yes).",/,
& "Default option values are:",/,
& "n = 1 (interpolate Gauss values to nodes)",/,
- & "a = 0 (do not average Gauss values to element centroids)",/,
- & "no = 1 (output values at nodes)",/,
- & "co = 0 (do not output values at element centroids)"
+ & "nf = 1 (output values at nodes)",/,
+ & "cf = 0 (do not output values at element centroids)")
c
return
end
c
c
+ subroutine fill(arr,val,len)
+c
+c... routine to fill a double-precision array with a given value
+c
+ implicit none
+ integer len
+ double precision arr(len),val
+c
+ integer i
+c
+ do i=1,len
+ arr(i)=val
+ end do
+ return
+ end
+c
+c
+ subroutine getjac(x,xs,det,shj,nsd,nen,iel)
+c
+c... subroutine to compute the jacobian determinant given the element
+c coordinates and the shape functions in natural coordinates.
+c
+c shj(1,nen),sh(2,nen),sh(3,nen) = x,y,and z derivatives
+c of shape functions
+c xs(nsd,nsd) = jacobian matrix
+c det = determinant of jacobian matrix
+c x(nsd,nen) = local nodal coordinates
+c
+ implicit none
+c
+c... subroutine arguments
+c
+ integer nsd,nen,iel
+ double precision x(nsd,nen),xs(nsd,nsd),det,shj(nsd+1,nen)
+c
+c... i/o units
+c
+ integer kti,kto,km,kp,kg,ko
+ common/units/kti,kto,km,kp,kg,ko
+c
+c...calculate jacobian matrix for (x,y,z) to (r,s,t) transformation
+c
+ call dgemm("n","t",nsd,nsd,nen,1.0d0,x,nsd,shj,nsd+1,0.0d0,xs,nsd)
+c
+c...form determinant of jacobian matrix and check for error condition
+c
+ det=xs(1,1)*xs(2,2)*xs(3,3)+xs(1,2)*xs(2,3)*xs(3,1)+xs(1,3)
+ & *xs(2,1)*xs(3,2)-xs(1,3)*xs(2,2)*xs(3,1)-xs(1,2)*xs(2,1)
+ & *xs(3,3)-xs(1,1)*xs(2,3)*xs(3,2)
+c
+ if(det.le.0.0d0) then
+ write(kto,700) iel,det
+ stop
+ end if
+c
+ 700 format("getjac: element # ",i7,2x,1pe15.8)
+ return
+ end
+c
+c
subroutine lcoord(ien,x,xl,nsd,nen,numnp)
c
c... subroutine to localize element coordinates
@@ -624,94 +559,235 @@
end
c
c
- subroutine pskip(iunit)
+ subroutine nodout(
+ & x,ien,numnp,numel,nsd,nen,ngauss,ninterp,
+ & vnodes,isnattr,nnattr,ntnattr,
+ & velemt,iseattr,neattr,nteattr,
+ & vnout,isnattrout,nnattrout,ntnattrout,nvnamesout,
+ & smop,sh,gauss,xl)
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... subroutine to interpolate Gauss point values to nodes
+c (if requested) and output values at nodes.
c
+ implicit none
c
+c... subroutine arguments
+c
+ integer numnp,numel,nsd,nen,ngauss
+ integer nnattr,ntnattr
+ integer neattr,nteattr
+ integer nnattrout,ntnattrout
+ integer ien(nen,numel)
+ integer isnattr(nnattr),iseattr(neattr),isnattrout(nnattrout)
+ double precision x(nsd,numnp)
+ double precision vnodes(ntnattr,numnp)
+ double precision velemt(nteattr,ngauss,numel)
+ double precision vnout(ntnattrout,numnp)
+ double precision smop(numnp),sh(nsd+1,nen,ngauss)
+ double precision gauss(nsd+1,ngauss),xl(nsd,nen)
+ logical ninterp
+ character nvnamesout(nnattrout)*(*)
+c
+c... i/o units
+c
+ integer kti,kto,km,kp,kg,ko
+ common/units/kti,kto,km,kp,kg,ko
+c
+c... local variables
+c
+ integer i,j,k,l,node
+ double precision xs(3,3),det
+c
+ if(ninterp) then
+ if(nen.eq.4) then
+ call plintet(sh,gauss,nsd,nen,ngauss)
+ else
+ call plinhex(sh,gauss,nsd,nen,ngauss)
+ end if
+c
+c... loop over elements and compute smoothing operator and weighted
+c state variables at nodes.
+c
+ call fill(smop,0.0d0,numnp)
+ do i=1,numel
+c
+c... localize nodal coordinates for element
+c
+ call lcoord(ien(1,i),x,xl,nsd,nen,numnp)
+c
+c... loop over Gauss points
+c
+ do j=1,ngauss
+ call getjac(xl,xs,det,sh,nsd,nen,i)
+ do k=1,nen
+ node=ien(k,i)
+ smop(node)=smop(node)+det
+ do l=1,nteattr
+ vnout(l,node)=vnout(l,node)+velemt(l,j,i)*det
+ end do
+ end do
+ end do
+ end do
+c
+c... invert smoothing operator and compute smoothed state variables
+c
+ do i=1,numnp
+ smop(i)=1.0d0/smop(i)
+ do j=1,nteattr
+ vnout(j,i)=vnout(j,i)*smop(i)
+ end do
+ end do
+ end if
+c
+c... output nodal values header
+c
+ write(ko,"(31i5)") nnattrout,(isnattrout(i),i=1,nnattrout)
+ do i=1,nnattrout
+ write(ko,"(a30)") nvnamesout(i)
+ end do
+c
+c... write nodal values
+c
+ if(.not.ninterp) then
+ do i=1,numnp
+ write(ko,"(i7,6(2x,1pe15.8))") i,(vnodes(j,i),j=1,ntnattr)
+ end do
+ else
+ do i=1,numnp
+ write(ko,"(i7,30(2x,1pe15.8))") i,(vnodes(j,i),j=1,ntnattr),
+ & (vnout(j,i),j=1,nteattr)
+ end do
+ end if
+ return
+ end
+c
+c
+ subroutine plinhex(sh,gauss,nsd,nen,ngauss)
+c
+c... Subroutine to compute shape functions in natural coordinates,
+c integration points, and weights for a trilinear hexahedron.
+c
implicit none
c
c... subroutine arguments
c
- integer iunit
+ integer nsd,nen,ngauss
+ double precision sh(nsd+1,nen,ngauss)
+ double precision gauss(nsd+1,ngauss)
c
c... local constants
c
- character leader*1
- data leader/'#'/
+ double precision r(8),s(8),t(8)
+ data r/-1d0, 1d0, 1d0,-1d0,-1d0, 1d0, 1d0,-1d0/
+ data s/-1d0,-1d0, 1d0, 1d0,-1d0,-1d0, 1d0, 1d0/
+ data t/ 1d0, 1d0, 1d0, 1d0,-1d0,-1d0,-1d0,-1d0/
c
c... intrinsic functions
c
- intrinsic index
+ intrinsic sqrt
c
-c... external functions
+c... user-defined functions
c
- integer nchar,nnblnk
- external nchar,nnblnk
c
c... local variables
c
- integer inblnk
- character string*80
+ integer i,l
+ double precision rr,ss,tt,root3i,eighth,one
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... definitions
c
+ one=1.0d0
+ root3i=sqrt(1.0d0/3.0d0)
+ eighth=1.0d0/8.0d0
c
- function tetcmp(xl)
+c... Linear hex definition
c
-c... program to compute determinant and equivalent tetrahedral volume
-c for a given 4x4 matrix.
+ do l=1,ngauss
+ gauss(1,l)=r(l)*root3i
+ gauss(2,l)=s(l)*root3i
+ gauss(3,l)=t(l)*root3i
+ gauss(4,l)=one
+ end do
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)
+ do l=1,ngauss
+ do i=1,nen
+ rr=one+r(i)*gauss(1,l)
+ ss=one+s(i)*gauss(2,l)
+ tt=one+t(i)*gauss(3,l)
+ sh(4,i,l)=eighth*rr*ss*tt
+ sh(1,i,l)=eighth*r(i)*ss*tt
+ sh(2,i,l)=eighth*s(i)*rr*tt
+ sh(3,i,l)=eighth*t(i)*rr*ss
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
+c
return
end
c
c
- function det3(x)
+ subroutine plintet(sh,gauss,nsd,nen,ngauss)
c
-c... function to compute determinant of a 3x3 matrix
+c... Subroutine to compute shape functions in natural coordinates,
+c integration points, and weights for a linear tetrahedron.
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)
+c
+c... subroutine arguments
+c
+ integer nsd,nen,ngauss
+ double precision sh(nsd+1,nen,ngauss)
+ double precision gauss(nsd+1,ngauss)
+c
+c... local constants
+c
+ double precision r(4),s(4),t(4),u(4)
+ data u/ 1d0, 0d0, 0d0, 0d0/
+ data r/ 0d0, 1d0, 0d0, 0d0/
+ data s/ 0d0, 0d0, 1d0, 0d0/
+ data t/ 0d0, 0d0, 0d0, 1d0/
+c
+c... intrinsic functions
+c
+c
+c... user-defined functions
+c
+c
+c... local variables
+c
+ integer i
+ double precision rr,ss,tt,uu
+ double precision tetvol,fourth
+c
+c... definitions
+c
+ tetvol=1.0d0/6.0d0
+ fourth=0.25d0
+c
+c... Linear hex definition
+c One-point integration is used in all cases.
+c
+ gauss(1,1)=fourth
+ gauss(2,1)=fourth
+ gauss(3,1)=fourth
+ gauss(4,1)=tetvol
+c
+ do i=1,nen
+ rr=r(i)*gauss(1,1)
+ ss=s(i)*gauss(2,1)
+ tt=t(i)*gauss(3,1)
+ uu=u(i)*gauss(3,1)
+ sh(4,i,1)=rr+ss+tt+uu
+ sh(1,i,1)=r(i)-u(i)
+ sh(2,i,1)=s(i)-u(i)
+ sh(3,i,1)=t(i)-u(i)
+ end do
+c
return
end
+c
+c version
+c $Id: plintet.f,v 1.4 2005/03/22 04:45:55 willic3 Exp $
+c
+c Generated automatically by Fortran77Mill on Wed May 21 14:15:03 2003
+c
+c End of file
More information about the cig-commits
mailing list