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