xref: /phasta/phSolver/common/hierarchic.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32) 
1c------------------------------------------------------------------------
2c
3c This file contains functions for dealing with higher order shape
4c functions at the element level.
5c
6c Christian Whiting, Winter 1999
7c------------------------------------------------------------------------
8
9      subroutine getsgn(ien, sgn)
10c------------------------------------------------------------------------
11c     returns the matrix of mode signs used for negating higher order
12c     basis functions. Connectivity array is assumed to have negative
13c     signs on all modes to be negated.
14c------------------------------------------------------------------------
15      include "common.h"
16
17      dimension ien(npro,nshl),  sgn(npro,nshl)
18
19      do i=nenl+1,nshl
20         where ( ien(:,i) < 0 )
21            sgn(:,i) = -one
22         elsewhere
23            sgn(:,i) = one
24         endwhere
25      enddo
26
27      return
28      end
29
30      subroutine getshp(shp, shgl, sgn, shape, shdrv)
31c------------------------------------------------------------------------
32c     returns the matrix of element shape functions with the higher
33c     order modes correctly negated at the current quadrature point.
34c------------------------------------------------------------------------
35      include "common.h"
36
37      dimension shp(nshl,ngauss),   shgl(nsd,nshl,ngauss),
38     &          sgn(npro,nshl),     shape(npro,nshl),
39     &          shdrv(npro,nsd,nshl)
40
41
42      do i=1,nenl
43         shape(:,i) = shp(i,intp)
44         do j=1,3
45            shdrv(:,j,i) = shgl(j,i,intp)
46         enddo
47      enddo
48      if ( ipord > 1 ) then
49         do i=nenl+1,nshl
50            shape(:,i) = sgn(:,i) * shp(i,intp)
51            do j=1,3
52               shdrv(:,j,i) = shgl(j,i,intp)*sgn(:,i)
53            enddo
54         enddo
55      endif
56
57      return
58      end
59
60      subroutine getshpb(shp, shgl, sgn, shape, shdrv)
61c------------------------------------------------------------------------
62c     returns the matrix of element shape functions with the higher
63c     order modes correctly negated at the current quadrature point.
64c------------------------------------------------------------------------
65      include "common.h"
66
67      dimension shp(nshl,ngaussb),  shgl(nsd,nshl,ngaussb),
68     &          sgn(npro,nshl),     shape(npro,nshl),
69     &          shdrv(npro,nsd,nshl)
70
71
72      do i=1,nenl
73         shape(:,i) = shp(i,intp)
74         do j=1,3
75            shdrv(:,j,i) = shgl(j,i,intp)
76         enddo
77      enddo
78      if ( ipord > 1 ) then
79         do i=nenl+1,nshl
80            shape(:,i) = sgn(:,i) * shp(i,intp)
81            do j=1,3
82               shdrv(:,j,i) = shgl(j,i,intp)*sgn(:,i)
83            enddo
84         enddo
85      endif
86
87      return
88      end
89
90      subroutine getbnodes(lnode)
91c------------------------------------------------------------------------
92c     compute the higher order modes that lie on the boundary of the
93c     element.
94c------------------------------------------------------------------------
95      include "common.h"
96
97      dimension lnode(27)
98
99c
100c.... boundary triangle of tet element
101c
102      if (lcsyst .eq. 1) then
103         do n = 1, nenbl
104            lnode(n) = n
105         enddo
106         if( ipord>1 ) then
107            nem = ipord-1
108            do n=1,3*nem
109               lnode(nenbl+n) = nenbl+1+n
110            enddo
111         endif
112         if( ipord>2 ) then
113            nfm = (ipord-2)*(ipord-1)/2
114            do n=1,nfm
115               lnode(3+3*nem+n) = 4+6*nem+n
116            enddo
117         endif
118c
119c.....boundary quadrilateral for a hex element
120c
121      else if(lcsyst .eq. 2) then
122         do n = 1, nenbl
123            lnode(n) = n
124         enddo
125         if( ipord > 1) then
126            nem = ipord -1
127            do n=1,4*nem
128               lnode(nenbl+n) = 8+n
129            enddo
130         endif
131
132         if( ipord > 3) then
133            nfm = (ipord-2)*(ipord-3)/2
134            do n=1,nfm
135               lnode(4+4*nem+n) = 8+12*nem+n
136            enddo
137         endif
138c
139c
140c.... This renumbers the boundary nodes for a wedge element when the
141c     boundary element is a quad face.  From lnode = [1 2 3 4]
142c                                        to   lnode = [1 4 5 2]
143c
144      else if(lcsyst .eq. 3) then
145         do n = 1, nenbl
146            lnode(n) = n
147         enddo
148c
149c     Need to implement for cubic, this valid only for ipord=2
150c
151         if( ipord>1 ) then
152            nem = ipord-1
153            do n=1,3*nem
154               lnode(nenbl+n) = 6+n
155            enddo
156         endif
157c
158c     Boundary quad of wedge element
159c
160      else if(lcsyst .eq. 4) then
161         lnode(1) = 1
162         lnode(2) = 4
163         lnode(3) = 5
164         lnode(4) = 2
165c$$$c
166c$$$c     Need to implement for cubic, this valid only for ipord=2
167c$$$c
168c$$$         if( ipord > 1) then
169c$$$            lnode(5) = 9
170c$$$            lnode(6) = 15
171c$$$            lnode(7) = 12
172c$$$            lnode(8) = 13
173c$$$               nem = ipord -1
174c$$$               do n=1,4*nem
175c$$$                  lnode(nenbl+n) = 6+n
176c$$$               enddo
177c$$$         endif
178c
179c     Boundary quad of pyramid element
180c
181      else if(lcsyst .eq. 5) then
182         lnode(1) = 1
183         lnode(2) = 2
184         lnode(3) = 3
185         lnode(4) = 4
186c$$$  c
187c$$$  c     Need to implement for cubic, this valid only for ipord=2
188c$$$  c
189c$$$            if( ipord > 1) then
190c$$$               lnode(5) = 9
191c$$$               lnode(6) = 15
192c$$$               lnode(7) = 12
193c$$$               lnode(8) = 13
194c$$$               nem = ipord -1
195c$$$               do n=1,4*nem
196c$$$                  lnode(nenbl+n) = 6+n
197c$$$               enddo
198c$$$            endif
199c
200c     Boundary triangle of pyramid element
201c
202      else if(lcsyst .eq. 6) then
203         lnode(1) = 1
204         lnode(2) = 5
205         lnode(3) = 2
206c$$$c
207c$$$c     Need to implement for cubic, this valid only for ipord=2
208c$$$c
209c$$$            if( ipord > 1) then
210c$$$               lnode(5) = 9
211c$$$               lnode(6) = 15
212c$$$               lnode(7) = 12
213c$$$               lnode(8) = 13
214c$$$               nem = ipord -1
215c$$$               do n=1,4*nem
216c$$$                  lnode(nenbl+n) = 6+n
217c$$$               enddo
218c$$$            endif
219c
220c.... other element types need to be implemented
221c
222      else
223         write (*,*) 'Boundary element not implemented for lcyst='
224     &        ,lcsyst
225         stop
226      endif
227
228      return
229      end
230
231c-----------------------------------------------------------------------
232c
233c  Evaluate coefficient vector at its interpolation points
234c
235c-----------------------------------------------------------------------
236      subroutine evalAtInterp( ycoeff,  yvals,  x,   nvars, npts )
237
238      use     pointer_data
239      include "common.h"
240
241      integer nvars, npts, nHits(nshg)
242
243      real*8  ycoeff(nshg,ndof),   yvals(nshg,nvars),
244     &        shp(nshl,npts),      shgl(nsd,nshl,npts),
245     &        intpnt(3,npts),      x(numnp,nsd)
246
247      real*8, allocatable :: ycl(:,:,:)
248      real*8, allocatable :: xl(:,:,:)
249      real*8, allocatable :: yvl(:,:,:)
250      real*8, allocatable :: sgn(:,:)
251
252      yvals = zero
253c
254c.... generate the shape functions at the interpolation points
255c
256      call getIntPnts(intpnt,npts)
257      do i=1,npts
258         call shpTet(ipord,intpnt(:,i),shp(:,i),shgl(:,:,i))
259      enddo
260c
261c.... loop over element blocks
262c
263      nHits = 0
264      do iblk = 1, nelblk
265         iel    = lcblk(1,iblk)
266         lcsyst = lcblk(3,iblk)
267         nenl   = lcblk(5,iblk) ! no. of vertices per element
268         nshl   = lcblk(10,iblk)
269         ndofl  = lcblk(8,iblk)
270         npro   = lcblk(1,iblk+1) - iel
271
272         allocate ( ycl(npro,nshl,ndof ) )
273         allocate ( yvl(npro,nshl,nvars) )
274         allocate ( xl(npro,nenl,nsd   ) )
275         allocate ( sgn(npro,nshl)       )
276
277         call getsgn(mien(iblk)%p,sgn)
278
279         call localy( ycoeff, ycl, mien(iblk)%p, ndof,  'gather  ')
280         call localx( x,      xl,  mien(iblk)%p, nsd,   'gather  ')
281
282         call eval  ( xl,       ycl,      yvl,
283     &                shp,      shgl,     sgn,
284     &                nvars,    npts    )
285
286c
287c.... average coefficients since stresses may be discontinuous
288c
289         call localSum( yvals,    yvl,    mien(iblk)%p,
290     &                  nHits,    nVars)
291
292
293         deallocate ( ycl )
294         deallocate ( yvl )
295         deallocate ( sgn )
296         deallocate ( xl  )
297c
298      enddo
299
300c
301c.... average the global values
302c
303      do i = 1, nshg
304         do j = 1, nvars
305            yvals(i,j) = yvals(i,j)/nHits(i) !(real(nHits(i),8))
306         enddo
307      enddo
308
309      return
310      end
311
312c-----------------------------------------------------------------------
313c
314c  evaluate in element coordinate system
315c
316c-----------------------------------------------------------------------
317      subroutine eval( xl,      ycl,     yvl,
318     &                 shp,     shgl,    sgn,
319     &                 nvars,   npts )
320
321      include "common.h"
322
323      integer nvars
324c
325      real*8  ycl(npro,nshl,ndof),   yvl(npro,nshl,nvars),
326     &        sgn(npro,nshl),        shape(npro,nshl),
327     &        shdrv(npro,nsd,nshl),  shp(nshl,npts),
328     &        shgl(nsd,nshl,npts),   xl(npro,nenl,nsd),
329     &        shg(npro,nshl,nsd),    gradV(npro,nsd,nsd),
330     &        dxidx(npro,nsd,nsd),   tmp(npro), wtmp
331
332      yvl = zero
333c
334c.... loop over interpolation points
335c
336      do intp = 1, npts
337         call getshp(shp,          shgl,      sgn,
338     &               shape,        shdrv)
339
340c
341c.... pressure and velocity
342c
343         do i = 1, nshl
344            do j = 1, 4
345               yvl(:,intp,j) = yvl(:,intp,j) + shape(:,i) * ycl(:,i,j)
346            enddo
347         enddo
348c
349c.... viscous stress
350c
351         call e3metric( xl,         shdrv,      dxidx,
352     &                  shg,        tmp)
353
354         gradV = zero
355         do n = 1, nshl
356            gradV(:,1,1) = gradV(:,1,1) + shg(:,n,1) * ycl(:,n,2)
357            gradV(:,2,1) = gradV(:,2,1) + shg(:,n,1) * ycl(:,n,3)
358            gradV(:,3,1) = gradV(:,3,1) + shg(:,n,1) * ycl(:,n,4)
359c
360            gradV(:,1,2) = gradV(:,1,2) + shg(:,n,2) * ycl(:,n,2)
361            gradV(:,2,2) = gradV(:,2,2) + shg(:,n,2) * ycl(:,n,3)
362            gradV(:,3,2) = gradV(:,3,2) + shg(:,n,2) * ycl(:,n,4)
363c
364            gradV(:,1,3) = gradV(:,1,3) + shg(:,n,3) * ycl(:,n,2)
365            gradV(:,2,3) = gradV(:,2,3) + shg(:,n,3) * ycl(:,n,3)
366            gradV(:,3,3) = gradV(:,3,3) + shg(:,n,3) * ycl(:,n,4)
367         enddo
368
369         rmu = datmat(1,2,1)
370
371         yvl(:,intp,6 ) = two * rmu * gradV(:,1,1)
372         yvl(:,intp,7 ) = two * rmu * gradV(:,2,2)
373         yvl(:,intp,8 ) = two * rmu * gradV(:,3,3)
374
375         yvl(:,intp,9 ) = rmu * ( gradV(:,1,2) + gradV(:,2,1) )
376         yvl(:,intp,10) = rmu * ( gradV(:,1,3) + gradV(:,3,1) )
377         yvl(:,intp,11) = rmu * ( gradV(:,2,3) + gradV(:,3,2) )
378
379c
380c.... loop over interpolation points
381c
382      enddo
383
384      return
385      end
386
387
388
389