xref: /phasta/phSolver/common/genint.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32) 
1      subroutine genint
2c
3c----------------------------------------------------------------------
4c
5c This subroutine inputs the integration information.
6c
7c----------------------------------------------------------------------
8c
9      include "common.h"
10
11      real*8, allocatable :: tmpQpt (:,:), tmpQwt (:)
12      real*8, allocatable :: tmpQptb(:,:), tmpQwtb(:)
13
14c
15c.... compute the shape function parameters
16c
17
18
19c
20c.... get quadrature data for interior and boundary elements
21c
22c
23c  Tets
24c
25      if (nen.eq.4) then
26         nshape  = (ipord+1)*(ipord+2)*(ipord+3)/6
27         nshapeb = (ipord+1)*(ipord+2)/2
28      endif
29
30      select case (intg(1,1))
31      case (1)
32         nint(1) = 1
33      case (2)
34         nint(1) = 4
35      case (3)
36         nint(1) = 16
37      case (4)
38         nint(1) = 29
39      end select
40      select case (intg(2,1))
41      case (1)
42         nintb(1) = 1
43      case (2)
44         nintb(1) = 3
45      case (3)
46         nintb(1) = 6
47      case (4)
48         nintb(1) = 12
49      end select
50
51      allocate (tmpQpt (4,nint(1)))
52      allocate (tmpQwt (nint(1)))
53      allocate (tmpQptb(4,nintb(1)))
54      allocate (tmpQwtb(nintb(1)))
55
56      call symtet(nint(1),tmpQpt,tmpQwt,nerr) ! interior elements
57      Qpt(1,1:4,1:nint(1)) = tmpQpt(1:4,1:nint(1))
58      Qwt(1,1:nint(1))     = tmpQwt(1:nint(1))
59
60
61      call symtri(nintb(1),tmpQptb,tmpQwtb,nerr) ! boundary elements
62      Qptb(1,1:4,1:nintb(1)) = tmpQptb(1:4,1:nintb(1))
63      Qwtb(1,1:nintb(1))     = tmpQwtb(1:nintb(1))
64
65      deallocate (tmpQpt)
66      deallocate (tmpQwt)
67      deallocate (tmpQptb)
68      deallocate (tmpQwtb)
69
70c
71c.... adjust quadrature weights to be consistent with the
72c     design of tau.
73c
74      Qwt(1,:) = (four/three)*Qwt(1,:)
75      Qwtb(1,:) = two*Qwtb(1,:)
76c
77c     Hexes now
78c
79      if (nen.eq.8) then
80         nshape = nen
81         if ( ipord .gt. 1 ) then
82            nshape = nshape + 12*(ipord - 1)
83         endif
84
85         if ( ipord .gt. 3 ) then
86            nshape = nshape + 3*(ipord -2)*(ipord - 3)
87         endif
88
89         if ( ipord .gt. 5 ) then
90            nshape = nshape + (ipord - 3)*(ipord - 4)*(ipord - 5)/6
91         endif
92
93         nshapeb = nenb
94
95         if ( ipord .gt. 1 ) then
96            nshapeb = nshapeb + 4*(ipord - 1)
97         endif
98
99         if ( ipord .gt. 3 ) then
100            nshapeb = nshapeb +(ipord -2)*(ipord - 3)/2
101         endif
102      endif
103
104      select case (intg(1,1))
105      case (1)
106         nint(2) = 1
107      case (2)
108         nint(2) = 8
109      case (3)
110         nint(2) =27
111      case (4)
112         nint(2) = 64
113      case (5)
114         nint(2) = 125
115      end select
116
117      select case (intg(2,1))
118      case (1)
119         nintb(2) = 1
120      case (2)
121         nintb(2) = 4
122      case(3)
123         nintb(2) = 9
124      case (4)
125         nintb(2) = 16
126      case (5)
127         nintb(2) = 25
128      end select
129
130
131      allocate (tmpQpt (4,nint(2)))
132      allocate (tmpQwt (nint(2)))
133      allocate (tmpQptb(4,nintb(2)))
134      allocate (tmpQwtb(nintb(2)))
135
136      call symhex(nint(2),tmpQpt,tmpQwt,nerr)
137      Qpt(2,1:4,1:nint(2)) = tmpQpt(1:4,1:nint(2))
138      Qwt(2,1:nint(2)) = tmpQwt(1:nint(2))
139
140      call symquad(nintb(2),tmpQptb,tmpQwtb,nerr)
141      Qptb(2,1:4,1:nintb(2)) = tmpQptb(1:4,1:nintb(2))
142      Qwtb(2,1:nintb(2)) = tmpQwtb(1:nintb(2))
143
144      deallocate (tmpQpt)
145      deallocate (tmpQwt)
146      deallocate (tmpQptb)
147      deallocate (tmpQwtb)
148
149      if (nen.eq.5) then        ! for pyramid
150
151         nshape = nen           ! counting the nodal functions
152
153         if ( ipord .gt. 1 ) then
154            nshape = nshape + 8*(ipord - 1) ! counting the edge functions
155         endif
156
157         if ( ipord .gt. 2 ) then
158                                ! counting the triangular face functions
159            nshape = nshape + 2*(ipord -1)*(ipord - 2)
160         endif
161
162         if ( ipord .gt. 3 ) then
163                                ! counting the quadrilateral face functions
164            nshape = nshape + (ipord - 2)*(ipord - 3)/2
165         endif
166
167         if ( ipord .gt. 5 ) then
168            nshape = nshape + (ipord - 3)*(ipord - 4)*(ipord - 5)/6
169         endif
170
171         nshapeb = nenb
172                                ! assume only the quadrilateral face on boundary
173         if ( ipord .gt. 1 ) then
174            nshapeb = nshapeb + 4*(ipord - 1)
175         endif
176         if ( ipord .gt. 3 ) then
177            nshape = nshape +(ipord - 2)*(ipord - 3)/2
178         endif
179
180      endif
181
182      select case (intg(1,1))
183      case (1)
184         nint(5) = 1
185      case (2)
186         nint(5) = 8
187      case (3)
188         nint(5) = 27
189      case (4)
190         nint(5) = 64
191      case (5)
192         nint(5) = 125
193      end select
194      select case (intg(2,1))
195      case (1)
196         nintb(5) = 1
197      case (2)
198         nintb(5) = 4
199      case (3)
200         nintb(5) = 9
201      case (4)
202         nintb(5) = 16
203      case (5)
204         nintb(5) = 25
205      end select
206      select case (intg(2,1))
207      case (1)
208         nintb(6) = 1
209      case (2)
210         nintb(6) = 4
211      case (3)
212         nintb(6) = 9
213      case (4)
214         nintb(6) = 12
215      case (5)
216         nintb(6) = 12 ! tri face boundary integration rules not coded
217      end select
218
219      allocate (tmpQpt (4,nint(5)))
220      allocate (tmpQwt (nint(5)))
221
222      call sympyr(nint(5),tmpQpt,tmpQwt,nerr) ! interior elements
223      Qpt(5,1:4,1:nint(5)) = tmpQpt(1:4,1:nint(5))
224      Qwt(5,1:nint(5)) = tmpQwt(1:nint(5))
225
226      deallocate (tmpQpt)
227      deallocate (tmpQwt)
228
229      allocate (tmpQptb(4,nintb(5)))
230      allocate (tmpQwtb(nintb(5)))
231
232      call symquad (nintb(5),tmpQptb,tmpQwtb,nerr) ! quad boundary elements
233      Qptb(5,1:4,1:nintb(5)) = tmpQptb(1:4,1:nintb(5))
234      Qwtb(5,1:nintb(5)) = tmpQwtb(1:nintb(5))
235
236      deallocate (tmpQptb)
237      deallocate (tmpQwtb)
238
239      allocate (tmpQptb(4,nintb(6)))
240      allocate (tmpQwtb(nintb(6)))
241
242      call symtripyr (nintb(6),tmpQptb,tmpQwtb,nerr) ! tri boundary elements
243      Qptb(6,1:4,1:nintb(6)) = tmpQptb(1:4,1:nintb(6))
244      Qwtb(6,1:nintb(6)) = tmpQwtb(1:nintb(6))
245
246      deallocate (tmpQptb)
247      deallocate (tmpQwtb)
248
249      if (nen.eq.6) then
250c
251c     later for hierarchic basis nshape formulas will be inserted
252c
253         nshape = nen
254
255         if ( ipord .gt. 1 ) then
256            nshape = nshape + 9*(ipord - 1)
257         endif
258
259         if ( ipord .gt. 2 ) then
260            nshape = nshape + (ipord -1)*(ipord - 2)
261         endif
262
263         if ( ipord .gt. 3 ) then
264            nshape = nshape + 3*(ipord - 2)*(ipord - 3)/2
265         endif
266
267         if ( ipord .gt. 4 ) then
268            nshape = nshape + (ipord - 2)*(ipord - 3)*(ipord - 4)/6
269         endif
270
271         nshapeb = nenb
272c
273         if (nenb .eq. 3) then
274            if ( ipord .gt. 1 ) then
275               nshapeb = nshapeb + 3*(ipord - 1)
276            endif
277            if ( ipord .gt. 2 ) then
278               nshape = nshape +(ipord - 1)*(ipord - 2)/2
279            endif
280         endif
281c
282         if (nenb .eq. 4) then
283            if ( ipord .gt. 1 ) then
284               nshapeb = nshapeb + 4*(ipord - 1)
285            endif
286            if ( ipord .gt. 2 ) then
287               nshape = nshape +(ipord - 2)*(ipord - 3)/2
288            endif
289         endif
290      endif
291
292      select case (intg(1,1))
293      case (2)
294         nint(3) = 6
295      case (3)
296         nint(3) = 18
297      case (4)
298         nint(3) = 48
299      end select
300      select case (intg(2,1))
301      case (2)
302         nintb(3) = 3
303         nintb(4) = 4
304      case (3)
305         nintb(3) = 6
306         nintb(4) = 9
307      case (4)
308         nintb(3) = 12
309         nintb(4) = 16
310      end select
311
312      allocate (tmpQpt (4,nint(3)))
313      allocate (tmpQwt (nint(3)))
314
315      call symwdg(nint(3),tmpQpt,tmpQwt,nerr) ! interior elements
316      Qpt(3,1:4,1:nint(3)) = tmpQpt(1:4,1:nint(3))
317      Qwt(3,1:nint(3)) = tmpQwt(1:nint(3))
318
319      deallocate (tmpQpt)
320      deallocate (tmpQwt)
321
322      allocate (tmpQptb(4,nintb(3)))
323      allocate (tmpQwtb(nintb(3)))
324
325      call symtri(nintb(3),tmpQptb,tmpQwtb,nerr) ! boundary elements
326      Qptb(3,1:2,2:nintb(3)) = tmpQptb(1:2,1:nintb(3)-1)
327      Qptb(3,1:2,1) = tmpQptb(1:2,nintb(3))
328c
329c     wedges want the third entry to be zeta=-1 (not t=1-r-s)
330c     4th entry not used
331c
332      Qptb(3,3:4,1:nintb(3)) = -1
333c$$$  Qptb(3,1:4,1:nintb(3)) = tmpQptb(1:4,1:nintb(3))
334      Qwtb(3,1:nintb(3)) = tmpQwtb(1:nintb(3))
335
336      deallocate (tmpQptb)
337      deallocate (tmpQwtb)
338
339      allocate (tmpQptb(4,nintb(4)))
340      allocate (tmpQwtb(nintb(4)))
341
342      call symquadw(nintb(4),tmpQptb,tmpQwtb,nerr) ! boundary elements
343      Qptb(4,1:4,1:nintb(4)) = tmpQptb(1:4,1:nintb(4))
344      Qwtb(4,1:nintb(4)) = tmpQwtb(1:nintb(4))
345
346      deallocate (tmpQptb)
347      deallocate (tmpQwtb)
348
349c
350c.... return
351c
352      return
353      end
354