Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ROUTINE Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine find_volumes (icon,octree,noctree,nleaves,lsf,nlsf,nnode,influid, &
x,y,z,u,v,w,s,n1,n2,n3,strain)
use gauss
implicit none
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! in s are the principal stresses (most compressive, least compressive, intermediary)
! and in n1,n2,n3 the corresponding directions
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
integer icon(8,nleaves)
integer octree(noctree)
integer noctree
integer nleaves
double precision lsf(nnode,nlsf)
integer nlsf
integer nnode
logical influid(nnode)
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer k,ie,ilsf,iint
integer,dimension(:),allocatable::levs
double precision phi(8),vol,voltot,work
double precision x(nnode),y(nnode),z(nnode),u(nnode),v(nnode),w(nnode)
double precision xx(8),yy(8),zz(8),vx(8),vy(8),vz(8)
double precision dhdx(8),dhdy(8),dhdz(8)
double precision volume,exx,eyy,ezz,exy,eyz,ezx
double precision s(3,nleaves),n1(3,nleaves),n2(3,nleaves),n3(3,nleaves)
double precision strain(3,3,nleaves)
logical todo
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
s=0.d0
n1=1.d-3
n2=1.d-3
n3=1.d-3
strain=0.d0
allocate (levs(nleaves))
call octree_find_element_level (octree,noctree,levs,nleaves)
do ilsf=1,nlsf
work=0.d0
voltot=0.d0
do ie=1,nleaves
exx=0.d0
eyy=0.d0
ezz=0.d0
exy=0.d0
eyz=0.d0
ezx=0.d0
do k=1,8
xx(k)=x(icon(k,ie))
yy(k)=y(icon(k,ie))
zz(k)=z(icon(k,ie))
vx(k)=u(icon(k,ie))
vy(k)=v(icon(k,ie))
vz(k)=w(icon(k,ie))
enddo
do iint=1,8
call compute_dhdx_dhdy_dhdz (8,rr(iint),ss(iint),tt(iint),xx,yy,zz,dhdx,dhdy,dhdz,volume)
do k=1,8
exx=exx+dhdx(k)*vx(k)
eyy=eyy+dhdy(k)*vy(k)
ezz=ezz+dhdz(k)*vz(k)
exy=exy+(dhdx(k)*vy(k)+dhdy(k)*vx(k))/2.d0
eyz=eyz+(dhdy(k)*vz(k)+dhdz(k)*vy(k))/2.d0
ezx=ezx+(dhdz(k)*vx(k)+dhdx(k)*vz(k))/2.d0
enddo
enddo
do k=1,8
phi(k)=lsf(icon(k,ie),ilsf)
enddo
call compute_positive_volume (phi,vol,3)
vol=1.d0-vol
vol=vol/(2.d0**levs(ie))**3
voltot=voltot+vol
work=work+(exx**2+eyy**2+ezz**2+2.d0*(exy**2+eyz**2+ezx**2))*vol
todo=.false.
do k=1,8
if (influid(icon(k,ie))) todo=.true.
enddo
!if (ilsf.eq.1.and.todo) &
if (todo) &
call PS (exx,eyy,ezz,exy,eyz,ezx,s(1,ie),s(2,ie),s(3,ie), &
n1(1,ie),n1(2,ie),n1(3,ie), &
n2(1,ie),n2(2,ie),n2(3,ie), &
n3(1,ie),n3(2,ie),n3(3,ie))
!if (ilsf.eq.1 .and. todo) then
if (todo) then
strain(1,1,ie)=exx
strain(1,2,ie)=exy
strain(1,3,ie)=ezx
strain(2,1,ie)=exy
strain(2,2,ie)=eyy
strain(2,3,ie)=eyz
strain(3,1,ie)=ezx
strain(3,2,ie)=eyz
strain(3,3,ie)=ezz
endif
enddo
write(*,'(a,i3)') 'ilsf=',ilsf
write(*,*) 'volume of material under: ',voltot
write(*,*) 'measured work: ',work
enddo
write(*,*) '--------------------------------------------------------------------------'
deallocate (levs)
return
end
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------