1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
*DECK SNQU01
SUBROUTINE SNQU01(NLF,JOP,U,W,TPQ,UPQ,VPQ,WPQ)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Set the level-symmetric (type 1) quadratures.
*
*Copyright:
* Copyright (C) 2005 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input
* NLF order of the SN approximation (even number).
*
*Parameters: output
* JOP number of base points per axial level in one octant.
* U base points in $\\xi$ of the axial quadrature. Used with
* zero-weight points.
* W weights for the axial quadrature in $\\xi$.
* TPQ base points in $\\xi$ of the 2D SN quadrature.
* UPQ base points in $\\mu$ of the 2D SN quadrature.
* VPQ base points in $\\eta$ of the 2D SN quadrature.
* WPQ weights of the 2D SN quadrature.
*
*-----------------------------------------------------------------------
*
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NLF,JOP(NLF/2)
REAL U(NLF/2),W(NLF/2),TPQ(NLF*(NLF/2+1)/4),UPQ(NLF*(NLF/2+1)/4),
1 VPQ(NLF*(NLF/2+1)/4),WPQ(NLF*(NLF/2+1)/4)
*----
* LOCAL VARIABLES
*----
PARAMETER(PI=3.141592654,MAXNLF=20,MAXNBA=55,MAXW=12)
INTEGER INWEI(MAXNBA)
DOUBLE PRECISION ZMU1,ZMU2,WSUM2,WEI(MAXW)
*----
* SET THE UNIQUE QUADRATURE VALUES.
*----
IF(NLF.GT.MAXNLF) CALL XABORT('SNQU01: MAXNLF OVERFLOW.')
M2=NLF/2
NPQ=M2*(M2+1)/2
NW=1+(NLF*(NLF+8)-1)/48
IF(NW.GT.MAXW) CALL XABORT('SNQU01: MAXW OVERFLOW.')
ZMU1=0.0D0
WEI(:NW)=0.0D0
IF(NLF.EQ.2) THEN
ZMU1=1.0D0/3.0D0
WEI(1)=1.0D0
ELSE IF(NLF.EQ.4) THEN
ZMU1=0.3500212**2
WEI(1)=1.0D0/3.0D0
ELSE IF(NLF.EQ.6) THEN
ZMU1=0.2666355**2
WEI(1)=0.1761263
WEI(2)=0.1572071
ELSE IF(NLF.EQ.8) THEN
ZMU1=1.0D0/21.0D0
WEI(1)=0.1209877
WEI(2)=0.0907407
WEI(3)=0.0925926
ELSE IF(NLF.EQ.12) THEN
ZMU1=0.1672126**2
WEI(1)=0.0707626
WEI(2)=0.0558811
WEI(3)=0.0373377
WEI(4)=0.0502819
WEI(5)=0.0258513
ELSE IF(NLF.EQ.16) THEN
ZMU1=0.1389568**2
WEI(1)=0.0489872
WEI(2)=0.0413296
WEI(3)=0.0212326
WEI(4)=0.0256207
WEI(5)=0.0360486
WEI(6)=0.0144586
WEI(7)=0.0344958
WEI(8)=0.0085179
ELSE
CALL XABORT('SNQU01: ORDER NOT AVAILABLE.')
ENDIF
U(1)=REAL(SQRT(ZMU1))
DO I=2,M2
ZMU2=ZMU1+2.0D0*DBLE(I-1)*(1.0D0-3.0D0*ZMU1)/DBLE(NLF-2)
U(I)=REAL(SQRT(ZMU2))
ENDDO
*----
* COMPUTE THE POSITION OF WEIGHTS.
*----
IPR=0
INMAX=0
DO IP=1,M2
JOP(IP)=M2-IP+1
DO IQ=1,JOP(IP)
IPR=IPR+1
IF(IPR.GT.MAXNBA) CALL XABORT('SNQU01: MAXNBA OVERFLOW.')
TPQ(IPR)=U(IP)
UPQ(IPR)=U(M2+2-IP-IQ)
VPQ(IPR)=U(IQ)
IS=MIN(IP,IQ,M2+2-IP-IQ)
NW0=0
DO II=1,IS-1
NW0=NW0+(M2-3*(II-1)+1)/2
ENDDO
KK=IP-IS+1
LL=IQ-IS+1
IF(KK.EQ.1)THEN
INWEI(IPR)=NW0+MIN(LL,M2-3*(IS-1)+1-LL)
ELSEIF(LL.EQ.1)THEN
INWEI(IPR)=NW0+MIN(KK,M2-3*(IS-1)+1-KK)
ELSE
INWEI(IPR)=NW0+MIN(KK,LL)
ENDIF
INMAX=MAX(INMAX,INWEI(IPR))
ENDDO
ENDDO
IF(INMAX.NE.NW) CALL XABORT('SNQU01: INVALID VALUE OD NW.')
IF(IPR.NE.NPQ) CALL XABORT('SNQU01: BAD VALUE ON NPQ.')
*----
* SET THE LEVEL-SYMMETRIC QUADRATURES.
*----
IPQ=0
WSUM=0.0
DO IP=1,M2
WSUM2=0.0D0
DO IQ=1,JOP(IP)
IPQ=IPQ+1
WPQ(IPQ)=REAL(WEI(INWEI(IPQ))*PI/2.0)
WSUM2=WSUM2+WEI(INWEI(IPQ))
ENDDO
W(IP)=REAL(WSUM2)
WSUM=WSUM+REAL(WSUM2*PI/2.0)
ENDDO
RETURN
END
|
