Newer
Older
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! MAKE_CUT Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
recursive subroutine make_cut (level,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,&
u,v,w,temp,pressure,strain,is_plastic,nnode,f, &
lsf,nlsf,r0,s0,t0,rst,icut,ileaves,eviscosity, &
vbounded,params,threadinfo,weightel)
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this subroutine is a intermediary routine between build_system and make_matrix
! to take into account the complex geometry of cut cells
! if we are in a non cut cell, make-matrix is called
! if we are in a cut cell but at a level that it smaller than levelmax, the
! cell is further cut and make_cut is recursively called
! if we are in a cut cell and level is equal to levelmax, we call make_matrix
! with material properties that have been interpolated from the
Dave Whipp
committed
! various material properties contribnuting to the cut cell
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
! level : current level in cut cell algorithm. It varies between 0 and levelmax.
! levelapprox : used to improve the postitive volutme calculation by further division
! mpe : number of nodes per element (8)
! ndof : number of degrees of freedom per node (3)
! ael : computed finite element matrix
! bel : computed rhs vector
! icon : connectivity array for the current element
! xg,yg,zg : global coordinate arrays of length nnode
! penalty : penalty factor used to impose the boundary conditions
! tempscale : temperature scaling parameter
! kfix : bc array of length ndof*nnode (kfix=1 means the dof is fixed to the value stored in velo)
! mat : material matrix for the nmat materials
! materialn : contains the material number associated to each lsf
! dt : time step length (only needed for the temperature calculations)
! u,v,w : velocity array (obtained from previous time step or at least containing the proper velocity at the fixed dofs)
! nnode : number of nodes
! f : global rhs vector
! lsf : global array of level set functions defining the surfaces
! nlsf : number of lsfs
! r0,s0,t0 : bottom left back corner of the part of the element. (we are now integrating in local r,s,t coordinates)
! rst : size of the part of the element we are integrating
! icut : returned, 0 if homogeneous element - 1 if cut element
! ileaves : current leaf number (useless except for debugging)
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use threads
use definitions
implicit none
type(parameters) params
integer level
integer levelmax
integer ndof
double precision ael(params%mpe*ndof,params%mpe*ndof)
double precision bel(params%mpe*ndof)
integer icon(params%mpe)
double precision x(nnode),y(nnode),z(nnode)
integer kfix(nnode*ndof)
type (material) mat(0:params%nmat)
double precision u(nnode),v(nnode),w(nnode)
double precision temp(nnode)
double precision pressure
double precision strain(nnode)
logical is_plastic
integer nnode
double precision f(nnode*ndof)
double precision lsf(params%mpe,nlsf)
integer nlsf
double precision r0,s0,t0,rst
integer icut
integer ileaves
double precision eviscosity
logical vbounded
type (thread) threadinfo
double precision weightel
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision volmax
double precision r(params%mpe),s(params%mpe),t(params%mpe)
double precision aelp(params%mpe*ndof,params%mpe*ndof),belp(params%mpe*ndof)
double precision h(params%mpe),vol_lsf0,prod
Dave Whipp
committed
double precision r0p,s0p,t0p,rstp,eps
double precision viscosity,density,penal,expon,diffusivity,heat,activ,expan
character (len=8) plasticity_type
double precision plasticity_parameters(9)
double precision,dimension(:,:),allocatable::lsfp
double precision,dimension(:),allocatable::vol_lsf
integer i,j,k,ii,jj,kk,jcut,levelp, err
integer matel
double precision eviscosityp
logical is_plastic_temp,vbounded_temp
double precision weight
Dave Whipp
committed
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
Dave Whipp
committed
eps=1.d-10
if (level.eq.0) then
ael=0.d0
bel=0.d0
eviscosity=0.d0
if (ndof.eq.3) weightel=0.d0
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
156
157
158
159
160
endif
matel=params%materialn(0)
if (.not.params%excl_vol) then
do i=1,nlsf
prod=lsf(1,i)
do k=2,params%mpe
if (prod*lsf(k,i).le.0.d0) goto 222
enddo
if (prod.lt.0.d0) then
matel=params%materialn(i)
endif
enddo
else
!check whether this is a cut cell
do i=1,nlsf
prod=lsf(1,i)
do k=2,params%mpe
if (prod*lsf(k,i).le.0.d0) goto 222
enddo
end do
!assign material to plain cell, since at that point we know that this is a plain cell
do i=1,nlsf
if (lsf(1,i).lt.0.d0) matel=params%materialn(i)
end do
!if (lsf(1,1).lt.0) matel=params%materialn(1)
!if (lsf(1,2).lt.0) matel=params%materialn(2)
end if
!=====[end new stuff]=====
Dave Whipp
committed
call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,mat(matel)%viscosity, &
mat(matel)%density,mat(matel)%penalty,mat(matel)%expon, &
mat(matel)%activationenergy,mat(matel)%expansion, &
mat(matel)%diffusivity,mat(matel)%heat, &
mat(matel)%plasticity_type,mat(matel)%plasticity_parameters, &
u,v,w,temp,pressure,strain,is_plastic_temp,nnode,f,r0,s0,t0, &
rst,ileaves,eviscosityp,vbounded_temp,threadinfo,weight)
ael=ael+aelp/(8.d0**level)
bel=bel+belp/(8.d0**level)
eviscosity=eviscosity+eviscosityp/(8.d0**level)
is_plastic=(is_plastic.or.is_plastic_temp)
vbounded=(vbounded.or.vbounded_temp)
if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
icut=0
return
!--------------------------------------------------------
222 continue
icut=1
! when we get to the bottom of the division we use an approximate algorithm
if (level.eq.levelmax) then
allocate (vol_lsf(nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc vol_lsf in make_cut$')
do i=1,nlsf
call compute_positive_volume (lsf(1,i),vol_lsf(i),params%levelapprox)
enddo
vol_lsf=1.d0-vol_lsf
if (.not.params%excl_vol) then
do i=1,nlsf-1
vol_lsf(i)=vol_lsf(i)-vol_lsf(i+1)
enddo
end if
! this is a little fix that is necessary due to the finite precision with
! which we compute the positive volumes and could lead to a contributing volume being
! either marginally negative or greater than 1
do i=1,nlsf
vol_lsf(i)=max(vol_lsf(i),0.d0)
vol_lsf(i)=min(vol_lsf(i),1.d0)
enddo
vol_lsf0=1.d0-sum(vol_lsf)
Dave Whipp
committed
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
if (vol_lsf0.gt.eps) then ! Element contains void, always use divFEM
viscosity=0.d0
density=0.d0
penal=0.d0
expon=0.d0
activ=0.d0
expan=0.d0
diffusivity=0.d0
heat=0.d0
do i=1,nlsf
matel=params%materialn(i)
viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
density=density+vol_lsf(i)*mat(matel)%density
penal=penal+vol_lsf(i)*mat(matel)%penalty
expon=expon+vol_lsf(i)*mat(matel)%expon
activ=activ+vol_lsf(i)*mat(matel)%activationenergy
expan=expan+vol_lsf(i)*mat(matel)%expansion
diffusivity=diffusivity+vol_lsf(i)/mat(matel)%diffusivity
heat=heat+vol_lsf(i)*mat(matel)%heat
enddo
matel=params%materialn(0)
viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
density=density+vol_lsf0*mat(matel)%density
penal=penal+vol_lsf0*mat(matel)%penalty
expon=expon+vol_lsf0*mat(matel)%expon
activ=activ+vol_lsf0*mat(matel)%activationenergy
expan=expan+vol_lsf0*mat(matel)%expansion
heat=heat+vol_lsf0*mat(matel)%heat
diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
diffusivity=1.d0/diffusivity ! Note that some properties add geometrically not algebraically
elseif (params%matrule.eq.1) then ! Assign material properties of volumetric majority material
matel=params%materialn(sum(maxloc(vol_lsf)))
viscosity=mat(matel)%viscosity
density=mat(matel)%density
penal=mat(matel)%penalty
expon=mat(matel)%expon
activ=mat(matel)%activationenergy
expan=mat(matel)%expansion
diffusivity=mat(matel)%diffusivity
heat=mat(matel)%heat
elseif (params%matrule.eq.2) then ! Assign material properties of volumetric minority material
matel=params%materialn(sum(minloc(vol_lsf)))
viscosity=mat(matel)%viscosity
density=mat(matel)%density
penal=mat(matel)%penalty
expon=mat(matel)%expon
activ=mat(matel)%activationenergy
expan=mat(matel)%expansion
diffusivity=mat(matel)%diffusivity
heat=mat(matel)%heat
else
if (params%matrule.ne.0) print *,'Invalid matrule value, using divFEM'
viscosity=0.d0
density=0.d0
penal=0.d0
expon=0.d0
activ=0.d0
expan=0.d0
diffusivity=0.d0
heat=0.d0
do i=1,nlsf
matel=params%materialn(i)
viscosity=viscosity+vol_lsf(i)*mat(matel)%viscosity
density=density+vol_lsf(i)*mat(matel)%density
penal=penal+vol_lsf(i)*mat(matel)%penalty
expon=expon+vol_lsf(i)*mat(matel)%expon
activ=activ+vol_lsf(i)*mat(matel)%activationenergy
expan=expan+vol_lsf(i)*mat(matel)%expansion
diffusivity=diffusivity+vol_lsf(i)/mat(matel)%diffusivity
heat=heat+vol_lsf(i)*mat(matel)%heat
enddo
matel=params%materialn(0)
viscosity=viscosity+vol_lsf0*mat(matel)%viscosity
density=density+vol_lsf0*mat(matel)%density
penal=penal+vol_lsf0*mat(matel)%penalty
expon=expon+vol_lsf0*mat(matel)%expon
activ=activ+vol_lsf0*mat(matel)%activationenergy
expan=expan+vol_lsf0*mat(matel)%expansion
heat=heat+vol_lsf0*mat(matel)%heat
diffusivity=diffusivity+vol_lsf0/mat(matel)%diffusivity
diffusivity=1.d0/diffusivity ! Note that some properties add geometrically not algebraically
endif
if (maxval(vol_lsf).gt.0.d0) then
if (params%matrule.eq.2) then
matel=params%materialn(sum(minloc(vol_lsf)))
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
else
matel=params%materialn(sum(maxloc(vol_lsf)))
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
endif
else
matel=params%materialn(0)
plasticity_type=mat(matel)%plasticity_type
plasticity_parameters=mat(matel)%plasticity_parameters
endif
!volmax=0.d0 ! Plasticity always uses a majority rule
!matel=params%materialn(0)
!plasticity_type=mat(matel)%plasticity_type
!plasticity_parameters=mat(matel)%plasticity_parameters
!do i=1,nlsf
! if (vol_lsf(i).gt.volmax) then
! volmax=vol_lsf(i)
! matel=params%materialn(i)
! plasticity_type=mat(matel)%plasticity_type
! plasticity_parameters=mat(matel)%plasticity_parameters
! endif
!enddo
! if (plasticity_type/='No') then
! print *,vol_lsf
! call stop_run ('pb$')
! end if
call make_matrix (params,ndof,aelp,belp,icon,x,y,z,kfix,viscosity,density, &
penal,expon,activ,expan,diffusivity,heat,plasticity_type,&
plasticity_parameters,u,v,w,temp,pressure,strain, &
is_plastic_temp,nnode,f,r0,s0,t0,rst,ileaves,eviscosityp,&
vbounded_temp,threadinfo,weight)
is_plastic=(is_plastic.or.is_plastic_temp)
vbounded=(vbounded.or.vbounded_temp)
ael=ael+aelp/(8.d0**level)
bel=bel+belp/(8.d0**level)
eviscosity=eviscosity+eviscosityp/(8.d0**level)
if (ndof.eq.3) weightel=weightel+weight/(8.d0**level)
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
deallocate (vol_lsf)
return
endif
!-----------------------------------------------------------------
! If we are not at the bottom level, we keep dividing
allocate (lsfp(params%mpe,nlsf),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lsfp in make_cut$')
do kk=1,2
t(1:4)=-1.d0+float(kk-1)
t(5:8)=float(kk-1)
do jj=1,2
s(1:2)=-1.d0+float(jj-1)
s(3:4)=float(jj-1)
s(5:6)=-1.d0+float(jj-1)
s(7:8)=float(jj-1)
do ii=1,2
r(1)=-1.d0+float(ii-1)
r(2)=float(ii-1)
r(3)=-1.d0+float(ii-1)
r(4)=float(ii-1)
r(5)=-1.d0+float(ii-1)
r(6)=float(ii-1)
r(7)=-1.d0+float(ii-1)
r(8)=float(ii-1)
do k=1,8
h(1)=(1.d0-r(k))*(1.d0-s(k))*(1.d0-t(k))/8.d0
h(2)=(1.d0+r(k))*(1.d0-s(k))*(1.d0-t(k))/8.d0
h(3)=(1.d0-r(k))*(1.d0+s(k))*(1.d0-t(k))/8.d0
h(4)=(1.d0+r(k))*(1.d0+s(k))*(1.d0-t(k))/8.d0
h(5)=(1.d0-r(k))*(1.d0-s(k))*(1.d0+t(k))/8.d0
h(6)=(1.d0+r(k))*(1.d0-s(k))*(1.d0+t(k))/8.d0
h(7)=(1.d0-r(k))*(1.d0+s(k))*(1.d0+t(k))/8.d0
h(8)=(1.d0+r(k))*(1.d0+s(k))*(1.d0+t(k))/8.d0
lsfp(k,:)=0.d0
do i=1,nlsf
do j=1,8
lsfp(k,i)=lsfp(k,i)+h(j)*lsf(j,i)
enddo
enddo
enddo
r0p=r0+rst*(r(1)+1.d0)/2.d0
s0p=s0+rst*(s(1)+1.d0)/2.d0
t0p=t0+rst*(t(1)+1.d0)/2.d0
rstp=rst/2.d0
levelp=level+1
Dave Whipp
committed
call make_cut (levelp,levelmax,ndof,ael,bel,icon,x,y,z,kfix,mat,u,v,w,&
temp,pressure,strain,is_plastic,nnode,f,lsfp,nlsf,r0p, &
s0p,t0p,rstp,jcut,ileaves,eviscosity,vbounded,params, &
threadinfo,weightel)
enddo
enddo
enddo
deallocate (lsfp)
return
end
!-------------------------------------------------------------------------------
Dave Whipp
committed
!-------------------------------------------------------------------------------