MFIX  2016-1
define_quadrics.f
Go to the documentation of this file.
1 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2 ! C
3 ! Module name: DEFINE_QUADRICS C
4 ! Purpose: Defines all matrices used to evaluate quadrics C
5 ! C
6 ! Author: Jeff Dietiker Date: 21-Feb-08 C
7 ! Reviewer: Date: C
8 ! C
9 ! Revision Number # Date: ##-###-## C
10 ! Author: # C
11 ! Purpose: # C
12 ! C
13 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
14  SUBROUTINE define_quadrics
15 
16  USE param
17  USE param1
18  USE parallel
19  USE constant
20  USE run
21  USE toleranc
22  USE geometry
23  USE indices
24  USE compar
25  USE sendrecv
26  USE quadric
27  USE cutcell
28  USE vtk
29 
30  IMPLICIT NONE
31  DOUBLE PRECISION, DIMENSION(3,3) :: Rx,Ry,Rz,C_QUADRIC,R_QUADRIC
32 
33 !======================================================================
34 ! Quadric Normal form : lambda_x * x^2 + lambda_y * y^2 + lambda_z * z^2 + d = 0
35 !======================================================================
36 
37  DO quadric_id = 1 , n_quadric
38 
39 ! Build translation matrices
40  CALL build_1x3_matrix(t_x(quadric_id),t_y(quadric_id),t_z(quadric_id),t_quadric(:,:,quadric_id))
41 ! Build characteristic matrices
42  CALL build_c_quadric_matrix(lambda_x(quadric_id),lambda_y(quadric_id),lambda_z(quadric_id),c_quadric)
43 ! Build Rotation matrices
44  CALL build_x_rotation_matrix(theta_x(quadric_id), rx)
45  CALL build_y_rotation_matrix(theta_y(quadric_id), ry)
46  CALL build_z_rotation_matrix(theta_z(quadric_id), rz)
47  r_quadric = matmul(rz,matmul(ry,rx))
48 ! Build A-matrices
49  a_quadric(:,:,quadric_id) = matmul(transpose(r_quadric),matmul(c_quadric,r_quadric))
50 
51  END DO
52 
53  ! Activate Quadric
54  quadric_id = n_quadric
55 
56  RETURN
57 
58  END SUBROUTINE define_quadrics
59 
60 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
61 ! C
62 ! Module name: GET_F_QUADRIC C
63 ! Purpose: Evaluate the function f(x,y,z) defining the quadric C
64 ! C
65 ! Author: Jeff Dietiker Date: 21-FEB-08 C
66 ! Reviewer: Date: C
67 ! C
68 ! C
69 ! Literature/Document References: C
70 ! C
71 ! Variables referenced: C
72 ! Variables modified: C
73 ! C
74 ! Local variables: C
75 ! C
76 ! Modified: ## Date: ##-###-## C
77 ! Purpose: ## C
78 ! C
79 ! Literature/Document References: ## C
80 ! C
81 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
82  SUBROUTINE get_f_quadric(x1,x2,x3,Q_ID,f,CLIP_FLAG)
83 
84  USE constant, only: pi
85  USE exit, only: mfix_exit
86  USE param1, only: half, one, undefined, zero
87  USE quadric
88  USE sendrecv, only: mype
89 
90  IMPLICIT NONE
91 
92  DOUBLE PRECISION :: x1,x2,x3,xt,yt,zt,R1,R2,xtr,ytr,ztr
93  DOUBLE PRECISION :: THETA,THETA1,THETA2,THETA3m,THETA3
94  DOUBLE PRECISION :: THETA1CYL1,THETA2TOR,THETA3CYL2,THETA_MIN
95  DOUBLE PRECISION :: Y1,Y2,R
96  DOUBLE PRECISION :: c,ss,ytest1,ytest2
97  DOUBLE PRECISION :: f
98  DOUBLE PRECISION :: fxmin,fxmax,fymin,fymax,fzmin,fzmax
99  DOUBLE PRECISION, DIMENSION(1,3) :: X_VECTOR,XMT
100  DOUBLE PRECISION, DIMENSION(3,1) :: TXMT
101  DOUBLE PRECISION, DIMENSION(1,1) :: TEMP_1x1
102  INTEGER :: Q_ID
103  LOGICAL :: CLIP_X,CLIP_Y,CLIP_Z,CLIP_FLAG
104  LOGICAL :: PIECE_X,PIECE_Y,PIECE_Z,PIECE_FLAG
105 
106  DOUBLE PRECISION :: YR1,YR2,RR1,RR2
107  DOUBLE PRECISION :: YRR1,YRR2,RC1,RC2,YC1,YC2
108 
109 
110 !======================================================================
111 
112  piece_x = (piece_xmin(q_id) <= x1).AND.( x1 <= piece_xmax(q_id))
113  piece_y = (piece_ymin(q_id) <= x2).AND.( x2 <= piece_ymax(q_id))
114  piece_z = (piece_zmin(q_id) <= x3).AND.( x3 <= piece_zmax(q_id))
115 
116  piece_flag = (piece_x.AND.piece_y.AND.piece_z)
117 
118  IF(.NOT.piece_flag) THEN
119  f = undefined
120  RETURN
121  ENDIF
122 
123  clip_x = (clip_xmin(q_id) <= x1).AND.( x1 <= clip_xmax(q_id))
124  clip_y = (clip_ymin(q_id) <= x2).AND.( x2 <= clip_ymax(q_id))
125  clip_z = (clip_zmin(q_id) <= x3).AND.( x3 <= clip_zmax(q_id))
126 
127  clip_flag = (clip_x.AND.clip_y.AND.clip_z)
128 
129 
130  IF(trim(quadric_form(q_id))=='PLANE') THEN
131 
132  f = lambda_x(q_id)*x1 + lambda_y(q_id)*x2 +lambda_z(q_id)*x3 + dquadric(q_id)
133 
134  ELSEIF(trim(quadric_form(q_id))=='TORUS_INT') THEN
135 
136  xt = x1-t_x(q_id)
137  yt = x2-t_y(q_id)
138  zt = x3-t_z(q_id)
139 
140  r1 = torus_r1(q_id)
141  r2 = torus_r2(q_id)
142 
143  f = -(4*(xt**2+zt**2)*r1**2-(xt**2+yt**2+zt**2+r1**2-r2**2)**2)
144 
145 
146  ELSEIF(trim(quadric_form(q_id))=='TORUS_EXT') THEN
147 
148  xt = x1-t_x(q_id)
149  yt = x2-t_y(q_id)
150  zt = x3-t_z(q_id)
151 
152  xtr = xt
153  ytr = yt
154  ztr = zt
155 
156  r1 = torus_r1(q_id)
157  r2 = torus_r2(q_id)
158 
159  f = 4*(xtr**2+ztr**2)*r1**2-(xtr**2+ytr**2+ztr**2+r1**2-r2**2)**2
160 
161 
162  ELSEIF(trim(quadric_form(q_id))=='Y_UCOIL_EXT') THEN
163 ! This shape represents a pair of parallel cylinders (y-direction)
164 ! capped at both ends by half a torus
165 ! to create a U-shaped coil
166 ! UCOIL_Y1 and UCOIL_Y2 are the min and max y-values of the coil
167 ! The coil is translated in the x and z direction (t_x, t_z)
168 ! and can be rotated about the y-direction, centered about (t_x,t_z)
169 ! by the angle THETA_Y
170 
171  c = dcos(theta_y(q_id))
172  ss = dsin(theta_y(q_id))
173 
174 ! Translation
175  xt = x1-t_x(q_id)
176  zt = x3-t_z(q_id)
177 
178 ! Rotation
179  xtr = xt*c + zt*ss
180  ztr = -xt*ss + zt*c
181 
182  r1 = ucoil_r1(q_id)
183  r2 = ucoil_r2(q_id)
184 
185 ! Limits of the cylinders.
186 ! There is half a torus above ytest1 and half a torus below ytest2
187  ytest1 = ucoil_y1(q_id) + r1 + r2
188  ytest2 = ucoil_y2(q_id) - (r1 + r2)
189 
190  IF(x2>=ytest2) THEN
191  ytr = x2 - ytest2
192  ELSEIF(x2<=ytest1) THEN
193  ytr = ytest1 - x2
194  ELSE
195  ytr = zero ! setting ytr = zero degenerates a torus into a pair of cylinders
196  ENDIF
197 
198  f = 4*(xtr**2+ytr**2)*r1**2-(xtr**2+ytr**2+ztr**2+r1**2-r2**2)**2
199 
200  ELSEIF(trim(quadric_form(q_id))=='Y_UCOIL2_EXT') THEN
201 ! This shape represents a pair of parallel cylinders (y-direction)
202 ! capped at both ends by a cylinder at 90 degree angle
203 ! to create a U-shaped coil
204 ! UCOIL_R1 is half the spacing between vertical cylinders
205 ! UCOIL_R2 is the cylinders radius
206 ! UCOIL_Y1 and UCOIL_Y2 are the min and max y-values of the coil
207 ! The coil is translated in the x and z direction (t_x, t_z)
208 ! and can be rotated about the y-direction, centered about (t_x,t_z)
209 ! by the angle THETA_Y
210 
211  c = dcos(theta_y(q_id))
212  ss = dsin(theta_y(q_id))
213 
214 ! Translation
215  xt = x1-t_x(q_id)
216  yt = x2
217  zt = x3-t_z(q_id)
218 
219 ! Rotation
220  xtr = xt*c + zt*ss
221  ytr = yt
222  ztr = -xt*ss + zt*c
223 
224  r1 = ucoil_r1(q_id)
225  r2 = ucoil_r2(q_id)
226 
227 
228 
229 
230 ! Limits of the cylinders.
231 ! There is half a torus above ytest1 and half a torus below ytest2
232  ytest1 = ucoil_y1(q_id) + r1 + r2
233  ytest2 = ucoil_y2(q_id) - (r1 + r2)
234 
235  IF(ytest1<=ytr.AND.ytr<=ytest2) THEN
236  f = 4*(xtr**2)*r1**2-(xtr**2+ztr**2+r1**2-r2**2)**2 ! a pair of cylinders
237  ELSEIF(ytr<=ytest1) THEN
238 ! Convert (x,y) into angle theta, and adjust its range from zero to 2*pi
239  theta = atan2(ytr-ytest1,xtr) ! Result is from -pi to pi
240  IF(-0.75*pi<=theta.AND.theta<=-0.25*pi) THEN
241  f = - ( ((ytr-(ucoil_y1(q_id) + r2))/r2)**2 + (ztr/r2)**2 -1.0 ) ! horizontal cylinder
242  ELSE
243  f = 4*(xtr**2)*r1**2-(xtr**2+ztr**2+r1**2-r2**2)**2 ! a pair of cylinders
244  ENDIF
245  ELSEIF(ytr>=ytest2) THEN
246 ! Convert (x,y) into angle theta, and adjust its range from zero to 2*pi
247  theta = atan2(ytr-ytest2,xtr) ! Result is from -pi to pi
248  IF(0.25*pi<=theta.AND.theta<0.75*pi) THEN
249  f = - ( ((ytr-(ucoil_y2(q_id) - r2))/r2)**2 + (ztr/r2)**2 -1.0 ) ! horizontal cylinder
250  ELSE
251  f = 4*(xtr**2)*r1**2-(xtr**2+ztr**2+r1**2-r2**2)**2 ! a pair of cylinders
252  ENDIF
253  ENDIF
254 
255 
256 
257  ELSEIF(trim(quadric_form(q_id))=='XY_BEND_INT') THEN
258 ! This shape represent a bend between two cylinders in the XY plane
259 ! BEND_R1 is the radius of the bend
260 ! BEND_R2 is the cylinders radius
261 ! BEND_THETA1 is the orientation of the first cylinder (Deg.)
262 ! BEND_THETA2 is the orientation of the second cylinder (Deg.).
263 ! The orientation is measured from the y-axis.
264 ! For example BEND_THETA1=0.0 and BEND_THETA2=90.0
265 ! represents a 90 deg bend between a vertical (first cylinder)
266 ! and a horizontal (second) cylinder.
267 ! The bend itself is represented by a section of a torus
268 ! The translation (t_x,t_y,t_z) defines the center of the bend (torus)
269 ! The shape is defines as three pieces: 2 cylinders and a torus
270 ! The switch between these pieces occur bases on angular position
271 ! There is a complicated definition of intermediate angles to allow
272 ! for a wide range of angles. This may be simplified in the future.
273 ! Both BEND_THETA1 and BEND_THETA2 must be between zero and 360.0 degrees.
274 ! Typically, BEND_THETA2 > BEND_THETA1, unless for a left bend,
275 ! For example, BEND_THETA1=330.0 and BEND_THETA2=30.0
276 ! will define a 60.0 deg. left-bend.
277 ! Specifying a bend angle larger that 180 degrees will likey fail
278 ! unless the bend is clipped somewhere else.
279 
280 ! Translation
281  xt = x1-t_x(q_id)
282  yt = x2-t_y(q_id)
283  zt = x3-t_z(q_id)
284 
285  r1 = bend_r1(q_id)
286  r2 = bend_r2(q_id)
287 ! Convert angles from degrees to radians
288  theta1 = bend_theta1(q_id)*(pi/180.0d0)
289  theta2 = bend_theta2(q_id)*(pi/180.0d0)
290 
291 ! Convert (x,y) into angle theta, and adjust its range from zero to 2*pi
292  theta = atan2(yt,xt) ! Result is from -pi to pi
293  IF(theta<zero) theta = theta + 2.0d0*pi
294 
295 ! THETA3 correspond to the point on a unit circle between THETA1 and THETA2
296  IF(theta2>theta1) THEN
297  theta3m = half*(theta1+theta2)
298  IF(theta3m<pi) THEN
299  theta3 = theta3m + pi
300  ELSE
301  theta3 = theta3m - pi
302  ENDIF
303  ELSE
304  theta3 = theta2 + half * (theta1-theta2)
305  ENDIF
306 
307 ! This angles are adjusted to wrap nicely along the discontinuity where
308 ! Theta=zero or 2*pi. There may be a simpler way to do this...
309 
310  IF(theta1>theta3) THEN
311  theta1cyl1 = theta1
312  ELSE
313  theta1cyl1 = theta1 + 2.0*pi
314  ENDIF
315 
316  IF(theta2>theta1) THEN
317  theta2tor = theta2
318  ELSE
319  theta2tor = theta2 + 2.0*pi
320  ENDIF
321 
322  IF(theta3>theta2) THEN
323  theta3cyl2 = theta3
324  ELSE
325  theta3cyl2 = theta3 + 2.0*pi
326  ENDIF
327 
328  theta_min = dmin1(theta1,theta2,theta3,theta1cyl1,theta2tor,theta3cyl2)
329 
330  IF(theta<theta_min) theta = theta + 2.0*pi
331 
332 ! Now join the pieces together:
333  IF(theta3<=theta.AND.theta<=theta1cyl1) THEN ! cylinder 1
334  c = dcos(theta1)
335  ss = dsin(theta1)
336 ! translation
337  xt = x1-(t_x(q_id)+r1*c)
338  yt = x2-(t_y(q_id)+r1*ss)
339  zt = x3-t_z(q_id)
340 ! Rotation
341  xtr = xt*c + yt*ss
342  ytr = -xt*ss + yt*c
343  ztr = zt
344 
345  f = (xtr/r2)**2 + (ztr/r2)**2 -1.0
346 
347  ELSEIF(theta1<=theta.AND.theta<=theta2tor) THEN ! Torus=elbow
348 ! translation
349  xt = x1-t_x(q_id)
350  yt = x2-t_y(q_id)
351  zt = x3-t_z(q_id)
352 ! Rotation
353  xtr = xt
354  ytr = yt
355  ztr = zt
356 
357  f = -4*(xtr**2+ytr**2)*r1**2+(xtr**2+ytr**2+ztr**2+r1**2-r2**2)**2
358 
359  ELSEIF(theta2<=theta.AND.theta<=theta3cyl2) THEN ! cylinder 2
360  c = dcos(theta2)
361  ss = dsin(theta2)
362 ! translation
363  xt = x1-(t_x(q_id)+r1*c)
364  yt = x2-(t_y(q_id)+r1*ss)
365  zt = x3-t_z(q_id)
366 ! Rotation
367  xtr = xt*c + yt*ss
368  ytr = -xt*ss + yt*c
369  ztr = zt
370 
371  f = (xtr/r2)**2 + (ztr/r2)**2 -1.0
372 
373  ELSE
374  WRITE(*,*)' Error in processing elbow.','THETA = ',theta/pi*180.0
375  CALL mfix_exit(mype)
376  ENDIF
377 
378  ELSEIF(trim(quadric_form(q_id))=='Y_C2C_INT') THEN
379 ! This shape connects two vertical cylinders by a conical section
380 ! This is a more convenient way of doing it the traditional way
381 ! with three quadrics (cylinder-cone-cylinder)
382 
383 ! Translation
384  xt = x1-t_x(q_id)
385  yt = x2-t_y(q_id)
386  zt = x3-t_z(q_id)
387 
388 ! Rotation
389  xtr = xt
390  ytr = yt
391  ztr = zt
392 
393 ! Radii
394  r1 = c2c_r1(q_id)
395  r2 = c2c_r2(q_id)
396 ! Extent of the conical section
397  y1 = c2c_y1(q_id)
398  y2 = c2c_y2(q_id)
399 
400  IF(ytr>=y2) THEN
401  r = r2
402  ELSEIF(ytr<=y1) THEN
403  r = r1
404  ELSE
405  IF(y2/=y1) THEN ! when Y2=Y1, then R2=R1 (see check_data_cartesian)
406  r = r1 + (r2-r1)/(y2-y1)*(ytr-y1)
407  ELSE
408  r = r1
409  ENDIF
410  ENDIF
411 
412  f = (xtr/r)**2 + (ztr/r)**2 - 1.0
413 
414 
415  ELSEIF(trim(quadric_form(q_id))=='REACTOR1') THEN
416 ! This shape defines a reactor (interior flow), made of two vertical cylinders,
417 ! connected by a conical section.
418 ! Each cylinder is rounded and closed by a conical cap.
419 
420 ! Translation
421  xt = x1-t_x(q_id)
422  yt = x2-t_y(q_id)
423  zt = x3-t_z(q_id)
424 
425 ! Rotation
426  xtr = xt
427  ytr = yt
428  ztr = zt
429 
430 ! Cylinders Radii
431  r1 = reactor1_r1(q_id) ! Lower section
432  r2 = reactor1_r2(q_id) ! Upper section
433 
434 ! Conical transition
435  y1 = reactor1_y1(q_id) ! Conical transition between lower
436  y2 = reactor1_y2(q_id) ! and upper sections
437 
438 ! Rounding
439  yr1 = reactor1_yr1(q_id) ! Lower rounding
440  rr1 = reactor1_rr1(q_id) ! Lower section rounding
441  theta1 = reactor1_theta1(q_id) ! angle (radians)
442 
443  yr2 = reactor1_yr2(q_id) ! Upper rounding
444  rr2 = reactor1_rr2(q_id) ! Lower section rounding
445  theta2 = reactor1_theta2(q_id) ! angle (radians)
446 
447 
448  yrr1 = yr1 - rr1 * dsin(theta1)
449  rc1 = r1-rr1*(one-dcos(theta1))
450  yc1 = yrr1 - rc1/dtan(theta1)
451  yrr2 = yr2 + rr2 * dsin(theta2)
452  rc2 = r2-rr2*(one-dcos(theta2))
453  yc2 = yrr2 + rc2/dtan(theta2)
454 
455  IF(ytr>=yc2) THEN ! Above upper conical cap
456 
457  r = zero
458 
459  ELSEIF(yrr2<=ytr.AND.ytr<=yc2) THEN ! upper conical cap
460 
461  r = rc2/(yrr2-yc2)*(ytr-yc2)
462 
463  ELSEIF(yr2<=ytr.AND.ytr<=yrr2) THEN ! upper rounding
464 
465  r = r2 - rr2 + dsqrt(rr2**2 - (yr2-ytr)**2)
466 
467  ELSEIF(y2<=ytr.AND.ytr<=yr2) THEN ! upper cylinder
468 
469  r = r2
470 
471  ELSEIF(y1<=ytr.AND.ytr<=y2) THEN ! conical transition
472 
473  r = r1 + (r2-r1)/(y2-y1)*(ytr-y1)
474 
475  ELSEIF(yr1<=ytr.AND.ytr<=y1) THEN ! lower cylinder
476 
477  r = r1
478 
479  ELSEIF(yrr1<=ytr.AND.ytr<=yr1) THEN ! lower rounding
480 
481  r = r1 - rr1 + dsqrt(rr1**2 - (yr1-ytr)**2)
482 
483  ELSEIF(yc1<=ytr.AND.ytr<=yrr1) THEN ! lower conical cap
484 
485  r = rc1/(yrr1-yc1)*(ytr-yc1)
486 
487  ELSE
488 
489  r = zero
490 
491  ENDIF
492 
493 
494  IF(r>zero) THEN
495  f = (xtr/r)**2 + (ztr/r)**2 - 1.0
496  ELSE
497  IF(ytr<=yc1) THEN
498  f = yc1 - ytr ! below lower conical cap
499  ELSE
500  f = ytr - yc2 ! above upper conical cap
501  ENDIF
502  ENDIF
503 
504 
505  ELSE
506 
507  CALL build_1x3_matrix(x1,x2,x3,x_vector)
508 
509  xmt = x_vector - t_quadric(:,:,q_id)
510  txmt = transpose(xmt)
511 
512  temp_1x1 = matmul(xmt,matmul(a_quadric(:,:,q_id),txmt))
513 
514  f = temp_1x1(1,1) + dquadric(q_id)
515 
516  ENDIF
517 
518 ! Each clipping limit is treated as a plane. For example, fxmin is
519 ! the equation of the plane describing x=xmin, and a value of fxmin
520 ! is compared with the current value of f to determine if the location
521 ! is part of the computational domain. The comparison (min of max)
522 ! follows the same logis as the 'AND' (max) , or 'OR' (min)
523 ! logic when combining two quadrics.
524 ! The clipping procedure is ignored when CLIP_FLAG is .FALSE.
525 ! This will happen when we are in a 'PIECEWISE' group
526 
527 
528 
529  IF(fluid_in_clipped_region(q_id)) THEN
530 
531  IF(clip_xmin(q_id)/=undefined) THEN
532  fxmin = -(clip_xmin(q_id)-x1)
533  f = dmin1(f,fxmin)
534  ENDIF
535 
536  IF(clip_xmax(q_id)/=undefined) THEN
537  fxmax = -(x1-clip_xmax(q_id))
538  f = dmin1(f,fxmax)
539  ENDIF
540 
541  IF(clip_ymin(q_id)/=undefined) THEN
542  fymin = -(clip_ymin(q_id)-x2)
543  f = dmin1(f,fymin)
544  ENDIF
545 
546  IF(clip_ymax(q_id)/=undefined) THEN
547  fymax = -(x2-clip_ymax(q_id))
548  f = dmin1(f,fymax)
549  ENDIF
550 
551  IF(clip_zmin(q_id)/=undefined) THEN
552  fzmin = -(clip_zmin(q_id)-x3)
553  f = dmin1(f,fzmin)
554  ENDIF
555 
556  IF(clip_zmax(q_id)/=undefined) THEN
557  fzmax = -(x3-clip_zmax(q_id))
558  f = dmin1(f,fzmax)
559  ENDIF
560 
561  ELSE
562 
563  IF(clip_xmin(q_id)/=undefined) THEN
564  fxmin = clip_xmin(q_id)-x1
565  f = dmax1(f,fxmin)
566  ENDIF
567 
568  IF(clip_xmax(q_id)/=undefined) THEN
569  fxmax = x1-clip_xmax(q_id)
570  f = dmax1(f,fxmax)
571  ENDIF
572 
573  IF(clip_ymin(q_id)/=undefined) THEN
574  fymin = clip_ymin(q_id)-x2
575  f = dmax1(f,fymin)
576  ENDIF
577 
578  IF(clip_ymax(q_id)/=undefined) THEN
579  fymax = x2-clip_ymax(q_id)
580  f = dmax1(f,fymax)
581  ENDIF
582 
583  IF(clip_zmin(q_id)/=undefined) THEN
584  fzmin = clip_zmin(q_id)-x3
585  f = dmax1(f,fzmin)
586  ENDIF
587 
588  IF(clip_zmax(q_id)/=undefined) THEN
589  fzmax = x3-clip_zmax(q_id)
590  f = dmax1(f,fzmax)
591  ENDIF
592 
593  ENDIF
594 
595  RETURN
596  END SUBROUTINE get_f_quadric
597 
598 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
599 ! C
600 ! Module name: REASSIGN_QUADRIC C
601 ! Purpose: Reassign the quadric based on location C
602 ! C
603 ! Author: Jeff Dietiker Date: 21-FEB-08 C
604 ! Reviewer: Date: C
605 ! C
606 ! C
607 ! Literature/Document References: C
608 ! C
609 ! Variables referenced: C
610 ! Variables modified: C
611 ! C
612 ! Local variables: C
613 ! C
614 ! Modified: ## Date: ##-###-## C
615 ! Purpose: ## C
616 ! C
617 ! Literature/Document References: ## C
618 ! C
619 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
620  SUBROUTINE reasssign_quadric(x1,x2,x3,GROUP,Q_ID)
621 
622  USE compar
623  USE exit, only: mfix_exit
624  USE quadric
625 
626  IMPLICIT NONE
627 
628  DOUBLE PRECISION x1,x2,x3
629  INTEGER :: I,Q_ID,GROUP,GS,P
630  LOGICAL :: PIECE_X,PIECE_Y,PIECE_Z,PIECE_FLAG
631  CHARACTER(LEN=9) :: GR
632 
633  q_id = 0
634 
635  gs = group_size(group)
636  gr = trim(group_relation(group))
637 
638  IF( gr /= 'PIECEWISE') RETURN
639 
640  DO p = 1 , gs
641 
642  i = group_q(group,p)
643 
644  piece_x = (piece_xmin(i) <= x1).AND.( x1 <= piece_xmax(i))
645  piece_y = (piece_ymin(i) <= x2).AND.( x2 <= piece_ymax(i))
646  piece_z = (piece_zmin(i) <= x3).AND.( x3 <= piece_zmax(i))
647 
648  piece_flag = (piece_x.AND.piece_y.AND.piece_z)
649 
650  IF (piece_flag) q_id = i
651 
652  ENDDO
653 
654  IF(q_id == 0 ) THEN
655  WRITE(*,*)' No Quadric defined at current location x,y,z=', x1,x2,x3
656  WRITE(*,*)' Please Check piecewise limits of quadric(s)'
657  WRITE(*,*)' Mfix will exit now.'
658  CALL mfix_exit(mype)
659  ENDIF
660 
661  RETURN
662 
663 
664  END SUBROUTINE reasssign_quadric
665 
666 
667 
668 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
669 ! C
670 ! Module name: BUILD_1x3_MATRIX C
671 ! Purpose: Catesian Grid - Build a (1x3) matrix C
672 ! from 3 scalars C
673 ! C
674 ! Author: Jeff Dietiker Date: 21-Feb-08 C
675 ! Reviewer: Date: C
676 ! C
677 ! Revision Number # Date: ##-###-## C
678 ! Author: # C
679 ! Purpose: # C
680 ! C
681 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
682  SUBROUTINE build_1x3_matrix(scalar1,scalar2,scalar3,M1x3)
683 
684  IMPLICIT NONE
685 
686  DOUBLE PRECISION:: scalar1,scalar2,scalar3
687  DOUBLE PRECISION, DIMENSION(1,3) :: M1x3
688  m1x3(1,1) = scalar1
689  m1x3(1,2) = scalar2
690  m1x3(1,3) = scalar3
691 
692  RETURN
693 
694  END SUBROUTINE build_1x3_matrix
695 
696 
697 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
698 ! C
699 ! Module name: BUILD_C_QUADRIC_MATRIX C
700 ! Purpose: Catesian Grid - Build a (3x3) diagonal matrix C
701 ! whose diagonal elements are the characteristic values of C
702 ! the quadric C
703 ! C
704 ! Author: Jeff Dietiker Date: 21-Feb-08 C
705 ! Reviewer: Date: C
706 ! C
707 ! Revision Number # Date: ##-###-## C
708 ! Author: # C
709 ! Purpose: # C
710 ! C
711 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
712  SUBROUTINE build_c_quadric_matrix(lambda1,lambda2,lambda3,C_QUADRIC)
713 
714  USE param1, only: zero
715 
716  IMPLICIT NONE
717 
718  DOUBLE PRECISION:: lambda1,lambda2,lambda3
719  DOUBLE PRECISION, DIMENSION(3,3) :: C_QUADRIC
720 
721 ! Transpose is used because matrices are stored column-wise
722 
723  c_quadric = transpose(reshape((/ &
724  lambda1 , zero , zero ,&
725  zero , lambda2 , zero ,&
726  zero , zero , lambda3 /),&
727  (/3,3/)))
728  RETURN
729 
730  END SUBROUTINE build_c_quadric_matrix
731 
732 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
733 ! C
734 ! Module name: BUILD_X_ROTATION_MATRIX C
735 ! Purpose: Catesian Grid - Build a (3x3) rotation matrix about x-axis,C
736 ! given the rotation angle in degrees C
737 ! C
738 ! Author: Jeff Dietiker Date: 21-Feb-08 C
739 ! Reviewer: Date: C
740 ! C
741 ! Revision Number # Date: ##-###-## C
742 ! Author: # C
743 ! Purpose: # C
744 ! C
745 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
746  SUBROUTINE build_x_rotation_matrix(Theta, R)
747 
748  USE param1, only: one, zero
749  USE constant, only: pi
750 
751  IMPLICIT NONE
752 
753  DOUBLE PRECISION:: Theta
754  DOUBLE PRECISION, DIMENSION(3,3) :: R
755 
756  theta = theta * (pi/180.0d0) ! Rotation angle about x-axis (radians)
757 
758 ! Transpose is used because matrices are stored column-wise
759 
760  r = transpose(reshape((/ &
761  one , zero , zero ,&
762  zero , dcos(theta) , dsin(theta) ,&
763  zero , -dsin(theta) , dcos(theta) /),&
764  (/3,3/)))
765 
766  RETURN
767 
768  END SUBROUTINE build_x_rotation_matrix
769 
770 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
771 ! C
772 ! Module name: BUILD_Y_ROTATION_MATRIX C
773 ! Purpose: Catesian Grid - Build a (3x3) rotation matrix about y-axis,C
774 ! given the rotation angle in degrees C
775 ! C
776 ! Author: Jeff Dietiker Date: 21-Feb-08 C
777 ! Reviewer: Date: C
778 ! C
779 ! Revision Number # Date: ##-###-## C
780 ! Author: # C
781 ! Purpose: # C
782 ! C
783 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
784  SUBROUTINE build_y_rotation_matrix(Theta, R)
785 
786  USE param1, only: one, zero
787  USE constant, only: pi
788 
789  IMPLICIT NONE
790 
791  DOUBLE PRECISION:: Theta
792  DOUBLE PRECISION, DIMENSION(3,3) :: R
793 
794  theta = theta * (pi/180.0d0) ! Rotation angle about x-axis (radians)
795 
796 ! Transpose is used because matrices are stored column-wise
797 
798  r = transpose(reshape((/ &
799  dcos(theta) , zero , dsin(theta) ,&
800  zero , one , zero ,&
801  -dsin(theta) , zero , dcos(theta) /),&
802  (/3,3/)))
803 
804  RETURN
805 
806  END SUBROUTINE build_y_rotation_matrix
807 
808 
809 !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
810 ! C
811 ! Module name: BUILD_Z_ROTATION_MATRIX C
812 ! Purpose: Catesian Grid - Build a (3x3) rotation matrix about z-axis,C
813 ! given the rotation angle in degrees C
814 ! C
815 ! Author: Jeff Dietiker Date: 21-Feb-08 C
816 ! Reviewer: Date: C
817 ! C
818 ! Revision Number # Date: ##-###-## C
819 ! Author: # C
820 ! Purpose: # C
821 ! C
822 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
823  SUBROUTINE build_z_rotation_matrix(Theta, R)
824 
825  USE param1, only: one, zero
826  USE constant, only: pi
827 
828  IMPLICIT NONE
829 
830  DOUBLE PRECISION:: Theta
831  DOUBLE PRECISION, DIMENSION(3,3) :: R
832 
833  theta = theta * (pi/180.0d0) ! Rotation angle about x-axis (radians)
834 
835 ! Transpose is used because matrices are stored column-wise
836 
837  r = transpose(reshape((/ &
838  dcos(theta) , dsin(theta) , zero ,&
839  -dsin(theta) , dcos(theta) , zero ,&
840  zero , zero , one /),&
841  (/3,3/)))
842 
843  RETURN
844 
845  END SUBROUTINE build_z_rotation_matrix
quadric::n_quadric
integer n_quadric
Definition: quadric_mod.f:6
quadric::t_quadric
double precision, dimension(1, 3, dim_quadric) t_quadric
Definition: quadric_mod.f:54
quadric::reactor1_y1
double precision, dimension(dim_quadric) reactor1_y1
Definition: quadric_mod.f:42
quadric::group_size
integer, dimension(dim_group) group_size
Definition: quadric_mod.f:68
param1::one
double precision, parameter one
Definition: param1_mod.f:29
quadric::theta_y
double precision, dimension(dim_quadric) theta_y
Definition: quadric_mod.f:20
quadric::c2c_r2
double precision, dimension(dim_quadric) c2c_r2
Definition: quadric_mod.f:36
quadric::c2c_r1
double precision, dimension(dim_quadric) c2c_r1
Definition: quadric_mod.f:36
quadric::lambda_x
double precision, dimension(dim_quadric) lambda_x
Definition: quadric_mod.f:14
quadric::torus_r2
double precision, dimension(dim_quadric) torus_r2
Definition: quadric_mod.f:24
get_f_quadric
subroutine get_f_quadric(x1, x2, x3, Q_ID, f, CLIP_FLAG)
Definition: define_quadrics.f:83
quadric::clip_xmax
double precision, dimension(dim_quadric) clip_xmax
Definition: quadric_mod.f:56
vtk
Definition: vtk_mod.f:1
quadric::t_x
double precision, dimension(dim_quadric) t_x
Definition: quadric_mod.f:18
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
quadric::ucoil_y1
double precision, dimension(dim_quadric) ucoil_y1
Definition: quadric_mod.f:28
quadric::ucoil_r2
double precision, dimension(dim_quadric) ucoil_r2
Definition: quadric_mod.f:26
quadric::reactor1_r1
double precision, dimension(dim_quadric) reactor1_r1
Definition: quadric_mod.f:40
quadric::c2c_y2
double precision, dimension(dim_quadric) c2c_y2
Definition: quadric_mod.f:34
quadric::torus_r1
double precision, dimension(dim_quadric) torus_r1
Definition: quadric_mod.f:24
exit::mfix_exit
subroutine mfix_exit(myID, normal_termination)
Definition: exit.f:5
quadric::theta_x
double precision, dimension(dim_quadric) theta_x
Definition: quadric_mod.f:20
exit
Definition: exit.f:2
quadric::ucoil_r1
double precision, dimension(dim_quadric) ucoil_r1
Definition: quadric_mod.f:26
quadric::bend_r2
double precision, dimension(dim_quadric) bend_r2
Definition: quadric_mod.f:30
quadric::bend_r1
double precision, dimension(dim_quadric) bend_r1
Definition: quadric_mod.f:30
quadric::reactor1_theta1
double precision, dimension(dim_quadric) reactor1_theta1
Definition: quadric_mod.f:48
quadric::dquadric
double precision, dimension(dim_quadric) dquadric
Definition: quadric_mod.f:16
param1::half
double precision, parameter half
Definition: param1_mod.f:28
run
Definition: run_mod.f:13
quadric::ucoil_y2
double precision, dimension(dim_quadric) ucoil_y2
Definition: quadric_mod.f:28
param
Definition: param_mod.f:2
quadric::lambda_y
double precision, dimension(dim_quadric) lambda_y
Definition: quadric_mod.f:14
define_quadrics
subroutine define_quadrics
Definition: define_quadrics.f:15
quadric::fluid_in_clipped_region
logical, dimension(dim_quadric) fluid_in_clipped_region
Definition: quadric_mod.f:60
quadric::quadric_id
integer quadric_id
Definition: quadric_mod.f:8
build_1x3_matrix
subroutine build_1x3_matrix(scalar1, scalar2, scalar3, M1x3)
Definition: define_quadrics.f:683
build_x_rotation_matrix
subroutine build_x_rotation_matrix(Theta, R)
Definition: define_quadrics.f:747
compar::mype
integer mype
Definition: compar_mod.f:24
quadric::quadric_form
character(len=12), dimension(dim_quadric) quadric_form
Definition: quadric_mod.f:10
build_z_rotation_matrix
subroutine build_z_rotation_matrix(Theta, R)
Definition: define_quadrics.f:824
quadric::group_relation
character(len=9), dimension(dim_group) group_relation
Definition: quadric_mod.f:72
build_y_rotation_matrix
subroutine build_y_rotation_matrix(Theta, R)
Definition: define_quadrics.f:785
quadric::clip_xmin
double precision, dimension(dim_quadric) clip_xmin
Definition: quadric_mod.f:56
build_c_quadric_matrix
subroutine build_c_quadric_matrix(lambda1, lambda2, lambda3, C_QUADRIC)
Definition: define_quadrics.f:713
quadric::reactor1_yr1
double precision, dimension(dim_quadric) reactor1_yr1
Definition: quadric_mod.f:44
quadric::t_y
double precision, dimension(dim_quadric) t_y
Definition: quadric_mod.f:18
quadric::piece_xmin
double precision, dimension(dim_quadric) piece_xmin
Definition: quadric_mod.f:58
quadric::reactor1_rr1
double precision, dimension(dim_quadric) reactor1_rr1
Definition: quadric_mod.f:46
constant::pi
double precision, parameter pi
Definition: constant_mod.f:158
quadric::bend_theta1
double precision, dimension(dim_quadric) bend_theta1
Definition: quadric_mod.f:32
quadric::t_z
double precision, dimension(dim_quadric) t_z
Definition: quadric_mod.f:18
quadric::a_quadric
double precision, dimension(3, 3, dim_quadric) a_quadric
Definition: quadric_mod.f:52
quadric::c2c_y1
double precision, dimension(dim_quadric) c2c_y1
Definition: quadric_mod.f:34
quadric::group_q
integer, dimension(dim_group, dim_quadric) group_q
Definition: quadric_mod.f:70
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
reasssign_quadric
subroutine reasssign_quadric(x1, x2, x3, GROUP, Q_ID)
Definition: define_quadrics.f:621
quadric::piece_xmax
double precision, dimension(dim_quadric) piece_xmax
Definition: quadric_mod.f:58