File: N:\mfix\model\cartesian_grid\define_quadrics.f

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
846