1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine printoutintfluidfem(prlab,ipkon,lakon,stx,
20     &  mi,ii,nelem,qfx,orab,ielorien,norien,co,kon,
21     &  ielmat,thicke,ielprop,prop,nelel,ithermal,orname,nelnew)
22!
23!     stores integration point results for element "nelem" in the .dat file
24!
25!     nelem is the element number from the input deck
26!     nelel is a local number, created e.g. by renumbering; thus
27!     far only applicable for FVM-CFD; for all other applications
28!     nelem=nelel
29!
30      implicit none
31!
32      character*6 prlab(*)
33      character*8 lakon(*)
34      character*80 orname(*)
35!
36      integer ipkon(*),mi(*),nelem,l,ii,mint3d,j,k,nope,
37     &  ielorien(mi(3),*),norien,kon(*),konl,indexe,m,iorien,iflag,
38     &  ielmat(mi(3),*),nopes,mint2d,kk,ki,kl,nlayer,ilayer,
39     &  null,ielprop(*),nelel,ithermal(*),nelef,nelnew(*)
40!
41      real*8 stx(6,mi(1),*),
42     &  qfx(3,mi(1),*),xi,et,ze,xl(3,20),xsj,shp(4,20),
43     &  coords(3,27),weight,orab(7,*),co(3,*),a(3,3),b(3,3),c(3,3),
44     &  qfxl(3),thicke(mi(3),*),xsj2(3),shp2(7,8),xl2(3,8),xs2(3,7),
45     &  thickness,tlayer(4),dlayer(4),xlayer(mi(3),4),
46     &  prop(*)
47!
48      include "gauss.f"
49!
50      data iflag /1/
51      data null /0/
52!
53      write(*,*) 'printoutintfluidfem ',nelel
54      write(*,*) 'printoutintfluidfem new ',nelnew(nelel)
55      nelef=nelnew(nelel)
56!
57      if(ipkon(nelef).lt.0) then
58         return
59      else
60         indexe=ipkon(nelef)
61      endif
62!
63!     check whether transformation is necessary (if orientation
64!     is applied and output in local system is requested)
65!
66      if(lakon(nelef)(7:8).ne.'LC') then
67         if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
68            iorien=0
69         else
70            iorien=max(0,ielorien(1,nelef))
71         endif
72      elseif(lakon(nelef)(4:5).eq.'20') then
73!
74!     composite materials
75!
76         if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
77            iorien=0
78         else
79!
80!           check whether at least one layer has a transformation
81!
82            iorien=0
83            do k=1,mi(3)
84               if(ielorien(k,nelef).ne.0) then
85                  iorien=max(0,ielorien(k,nelef))
86                  exit
87               endif
88            enddo
89         endif
90!
91!     determining the number of layers
92!
93         mint2d=4
94         nopes=8
95!
96         nlayer=0
97         do k=1,mi(3)
98            if(ielmat(k,nelef).ne.0) then
99               nlayer=nlayer+1
100            endif
101         enddo
102!
103!     determining the layer thickness and global thickness
104!     at the shell integration points
105!
106         iflag=1
107         do kk=1,mint2d
108            xi=gauss3d2(1,kk)
109            et=gauss3d2(2,kk)
110            call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
111            tlayer(kk)=0.d0
112            do k=1,nlayer
113               thickness=0.d0
114               do j=1,nopes
115                  thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
116               enddo
117               tlayer(kk)=tlayer(kk)+thickness
118               xlayer(k,kk)=thickness
119            enddo
120         enddo
121         iflag=3
122!
123         ilayer=0
124         do k=1,4
125            dlayer(k)=0.d0
126         enddo
127!
128!
129      elseif(lakon(nelef)(4:5).eq.'15') then
130!
131!     composite materials
132!
133         if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
134            iorien=0
135         else
136!
137!           check whether at least one layer has a transformation
138!
139            iorien=0
140            do k=1,mi(3)
141               if(ielorien(k,nelef).ne.0) then
142                  iorien=max(0,ielorien(k,nelef))
143                  exit
144               endif
145            enddo
146         endif
147!
148!     determining the number of layers
149!
150         mint2d=3
151         nopes=6
152!
153         nlayer=0
154         do k=1,mi(3)
155            if(ielmat(k,nelef).ne.0) then
156               nlayer=nlayer+1
157            endif
158         enddo
159!
160!     determining the layer thickness and global thickness
161!     at the shell integration points
162!
163         iflag=1
164         do kk=1,mint2d
165            xi=gauss3d10(1,kk)
166            et=gauss3d10(2,kk)
167            call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
168            tlayer(kk)=0.d0
169            do k=1,nlayer
170               thickness=0.d0
171               do j=1,nopes
172                  thickness=thickness+thicke(k,indexe+j)*shp2(4,j)
173               enddo
174               tlayer(kk)=tlayer(kk)+thickness
175               xlayer(k,kk)=thickness
176            enddo
177         enddo
178         iflag=3
179!
180         ilayer=0
181         do k=1,3
182            dlayer(k)=0.d0
183         enddo
184!
185      endif
186!
187!     number of integration points
188!
189      if(lakon(nelef)(4:5).eq.'8R') then
190         mint3d=1
191      elseif(lakon(nelef)(4:7).eq.'20RB') then
192         if((lakon(nelef)(8:8).eq.'R').or.
193     &      (lakon(nelef)(8:8).eq.'C')) then
194            mint3d=50
195         else
196            call beamintscheme(lakon(nelef),mint3d,ielprop(nelel),prop,
197     &           null,xi,et,ze,weight)
198         endif
199      elseif((lakon(nelef)(4:4).eq.'8').or.
200     &       (lakon(nelef)(4:6).eq.'20R')) then
201         if(lakon(nelef)(7:8).eq.'LC') then
202            mint3d=8*nlayer
203         else
204            mint3d=8
205         endif
206      elseif(lakon(nelef)(4:4).eq.'2') then
207         mint3d=27
208      elseif(lakon(nelef)(4:5).eq.'10') then
209         mint3d=4
210      elseif(lakon(nelef)(4:4).eq.'4') then
211         mint3d=1
212      elseif(lakon(nelef)(4:5).eq.'15') then
213         if(lakon(nelef)(7:8).eq.'LC') then
214            mint3d=6*nlayer
215         else
216            mint3d=9
217         endif
218      elseif(lakon(nelef)(4:4).eq.'6') then
219         mint3d=2
220      elseif(lakon(nelef)(1:1).eq.'U') then
221         mint3d=ichar(lakon(nelef)(6:6))
222      else
223         return
224      endif
225!
226!     calculation of the integration point coordinates for
227!     output in the local system (if needed)
228!
229      if((iorien.ne.0).or.(prlab(ii)(1:4).eq.'COOR')) then
230         if(lakon(nelef)(4:4).eq.'2') then
231            nope=20
232         elseif(lakon(nelef)(4:4).eq.'8') then
233            nope=8
234         elseif(lakon(nelef)(4:5).eq.'10') then
235            nope=10
236         elseif(lakon(nelef)(4:4).eq.'4') then
237            nope=4
238         elseif(lakon(nelef)(4:5).eq.'15') then
239            nope=15
240         elseif(lakon(nelef)(4:4).eq.'6') then
241            nope=6
242         endif
243!
244         do j=1,nope
245            konl=kon(indexe+j)
246            do k=1,3
247               xl(k,j)=co(k,konl)
248            enddo
249         enddo
250!
251         do j=1,mint3d
252            if((lakon(nelef)(4:5).eq.'8R').or.
253     &         (lakon(nelef)(1:4).eq.'F3D8')) then
254               xi=gauss3d1(1,j)
255               et=gauss3d1(2,j)
256               ze=gauss3d1(3,j)
257               weight=weight3d1(j)
258            elseif(lakon(nelef)(4:8).eq.'20RB') then
259               if((lakon(nelef)(8:8).eq.'B').or.
260     &            (lakon(nelef)(8:8).eq.'C')) then
261                  xi=gauss3d13(1,j)
262                  et=gauss3d13(2,j)
263                  ze=gauss3d13(3,j)
264                  weight=weight3d13(j)
265               else
266                  call beamintscheme(lakon(nelef),mint3d,ielprop(nelel),
267     &                  prop,j,xi,et,ze,weight)
268               endif
269            elseif((lakon(nelef)(4:4).eq.'8').or.
270     &             (lakon(nelef)(4:6).eq.'20R'))
271     &              then
272               if(lakon(nelef)(7:8).ne.'LC') then
273                  xi=gauss3d2(1,j)
274                  et=gauss3d2(2,j)
275                  ze=gauss3d2(3,j)
276                  weight=weight3d2(j)
277               else
278                  kl=mod(j,8)
279                  if(kl.eq.0) kl=8
280!
281                  xi=gauss3d2(1,kl)
282                  et=gauss3d2(2,kl)
283                  ze=gauss3d2(3,kl)
284                  weight=weight3d2(kl)
285!
286                  ki=mod(j,4)
287                  if(ki.eq.0) ki=4
288!
289                  if(kl.eq.1) then
290                     ilayer=ilayer+1
291                     if(ilayer.gt.1) then
292                        do k=1,4
293                           dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
294                        enddo
295                     endif
296                  endif
297                  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
298     &                 tlayer(ki)-1.d0
299                  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
300               endif
301            elseif(lakon(nelef)(4:4).eq.'2') then
302               xi=gauss3d3(1,j)
303               et=gauss3d3(2,j)
304               ze=gauss3d3(3,j)
305               weight=weight3d3(j)
306            elseif(lakon(nelef)(4:5).eq.'10') then
307               xi=gauss3d5(1,j)
308               et=gauss3d5(2,j)
309               ze=gauss3d5(3,j)
310               weight=weight3d5(j)
311            elseif(lakon(nelef)(4:4).eq.'4') then
312               xi=gauss3d4(1,j)
313               et=gauss3d4(2,j)
314               ze=gauss3d4(3,j)
315               weight=weight3d4(j)
316            elseif(lakon(nelef)(4:5).eq.'15') then
317               if(lakon(nelef)(7:8).ne.'LC') then
318                  xi=gauss3d8(1,j)
319                  et=gauss3d8(2,j)
320                  ze=gauss3d8(3,j)
321                  weight=weight3d8(j)
322               else
323                  kl=mod(j,6)
324                  if(kl.eq.0) kl=6
325!
326                  xi=gauss3d10(1,kl)
327                  et=gauss3d10(2,kl)
328                  ze=gauss3d10(3,kl)
329                  weight=weight3d10(kl)
330!
331                  ki=mod(j,3)
332                  if(ki.eq.0) ki=3
333!
334                  if(kl.eq.1) then
335                     ilayer=ilayer+1
336                     if(ilayer.gt.1) then
337                        do k=1,3
338                           dlayer(k)=dlayer(k)+xlayer(ilayer-1,k)
339                        enddo
340                     endif
341                  endif
342                  ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/
343     &                 tlayer(ki)-1.d0
344                  weight=weight*xlayer(ilayer,ki)/tlayer(ki)
345               endif
346            elseif(lakon(nelef)(1:4).eq.'C3D6') then
347               xi=gauss3d7(1,j)
348               et=gauss3d7(2,j)
349               ze=gauss3d7(3,j)
350               weight=weight3d7(j)
351            elseif(lakon(nelef)(1:4).eq.'F3D6') then
352               xi=gauss3d14(1,j)
353               et=gauss3d14(2,j)
354               ze=gauss3d14(3,j)
355               weight=weight3d14(j)
356            endif
357!
358            if(nope.eq.20) then
359               call shape20h(xi,et,ze,xl,xsj,shp,iflag)
360            elseif(nope.eq.8) then
361               call shape8h(xi,et,ze,xl,xsj,shp,iflag)
362            elseif(nope.eq.10) then
363               call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
364            elseif(nope.eq.4) then
365               call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
366            elseif(nope.eq.15) then
367               call shape15w(xi,et,ze,xl,xsj,shp,iflag)
368            else
369               call shape6w(xi,et,ze,xl,xsj,shp,iflag)
370            endif
371!
372            do k=1,3
373               coords(k,j)=0.d0
374               do l=1,nope
375                  coords(k,j)=coords(k,j)+xl(k,l)*shp(4,l)
376               enddo
377            enddo
378         enddo
379      endif
380!
381      if((prlab(ii)(1:4).eq.'S   ').or.(prlab(ii)(1:4).eq.'SVF ')) then
382         do j=1,mint3d
383!
384!           composite materials
385!
386            if(lakon(nelef)(7:8).eq.'LC') then
387               if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
388                  iorien=0
389               elseif(lakon(nelef)(4:5).eq.'20') then
390                  kl=mod(j,8)
391                  if(kl.eq.0) kl=8
392                  ilayer=(j-kl)/8+1
393                  iorien=max(0,ielorien(ilayer,nelef))
394               elseif(lakon(nelef)(4:5).eq.'15') then
395                  kl=mod(j,6)
396                  if(kl.eq.0) kl=6
397                  ilayer=(j-kl)/6+1
398                  iorien=max(0,ielorien(ilayer,nelef))
399               endif
400            endif
401!
402            if(iorien.eq.0) then
403               write(5,'(i10,1x,i3,1p,6(1x,e13.6))') nelem,j,
404     &              (stx(k,j,nelef),k=1,6)
405            else
406               call transformatrix(orab(1,iorien),coords(1,j),a)
407               b(1,1)=stx(1,j,nelef)
408               b(2,2)=stx(2,j,nelef)
409               b(3,3)=stx(3,j,nelef)
410               b(1,2)=stx(4,j,nelef)
411               b(1,3)=stx(5,j,nelef)
412               b(2,3)=stx(6,j,nelef)
413               b(2,1)=b(1,2)
414               b(3,1)=b(1,3)
415               b(3,2)=b(2,3)
416               do k=1,3
417                  do l=1,3
418                     c(k,l)=0.d0
419                     do m=1,3
420                        c(k,l)=c(k,l)+b(k,m)*a(m,l)
421                     enddo
422                  enddo
423               enddo
424               do k=1,3
425                  do l=k,3
426                     b(k,l)=0.d0
427                     do m=1,3
428                        b(k,l)=b(k,l)+a(m,k)*c(m,l)
429                     enddo
430                  enddo
431               enddo
432               write(5,'(i10,1x,i3,1p,6(1x,e13.6),1x,a20)') nelem,j,
433     &              b(1,1),b(2,2),b(3,3),b(1,2),b(1,3),b(2,3),
434     &              orname(iorien)(1:20)
435            endif
436         enddo
437      elseif(((prlab(ii)(1:4).eq.'HFL ').or.(prlab(ii)(1:4).eq.'HFLF'))
438     &        .and.(ithermal(1).gt.1)) then
439         do j=1,mint3d
440!
441!           composite materials
442!
443            if(lakon(nelef)(7:8).eq.'LC') then
444               if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then
445                  iorien=0
446               elseif(lakon(nelef)(4:5).eq.'20') then
447                  kl=mod(j,8)
448                  if(kl.eq.0) kl=8
449                  ilayer=(j-kl)/8+1
450                  iorien=max(0,ielorien(ilayer,nelef))
451               elseif(lakon(nelef)(4:5).eq.'15') then
452                  kl=mod(j,6)
453                  if(kl.eq.0) kl=6
454                  ilayer=(j-kl)/6+1
455                  iorien=max(0,ielorien(ilayer,nelef))
456               endif
457            endif
458!
459            if(iorien.eq.0) then
460               write(5,'(i10,1x,i3,1p,3(1x,e13.6))') nelem,j,
461     &              (qfx(k,j,nelef),k=1,3)
462            else
463               do k=1,3
464                  qfxl(k)=qfx(k,j,nelef)
465               enddo
466               call transformatrix(orab(1,iorien),coords(1,j),a)
467               write(5,'(i10,1x,i3,1p,3(1x,e13.6),1x,a20)') nelem,j,
468     &              qfxl(1)*a(1,1)+qfxl(2)*a(2,1)+qfxl(3)*a(3,1),
469     &              qfxl(1)*a(1,2)+qfxl(2)*a(2,2)+qfxl(3)*a(3,2),
470     &              qfxl(1)*a(1,3)+qfxl(2)*a(2,3)+qfxl(3)*a(3,3),
471     &              orname(iorien)(1:20)
472            endif
473         enddo
474      elseif(prlab(ii)(1:4).eq.'COOR') then
475         do j=1,mint3d
476            write(5,'(i10,1x,i3,1p,3(1x,e13.6))') nelem,j,
477     &           (coords(k,j),k=1,3)
478         enddo
479      endif
480!
481      return
482      end
483