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