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 calcexternalwork(co,vold,istartset,iendset,ipkon,
20     &  lakon,kon,ialset,nset,set,mi,externalfaces,nelemload,
21     &  nload,xload,sideload,delexternalwork,voldprev)
22!
23!     calculation and printout of the lift and drag forces
24!
25      implicit none
26!
27      character*8 lakonl,lakon(*)
28      character*20 sideload(*)
29      character*81 set(*),externalfaces
30!
31      integer konl(20),ifaceq(8,6),nelem,i,j,i1,j1,
32     &  jj,ig,nope,nopes,igpres,id,
33     &  mint2d,ifacet(6,4),ifacew(8,5),iflag,indexe,jface,istartset(*),
34     &  iendset(*),ipkon(*),kon(*),iset,ialset(*),nset,ipos,
35     &  mi(*),nelemload(2,*),nload,ipres
36!
37      real*8 co(3,*),xs2(3,7),pres,voldl2(3,8),xl2(3,8),xsj2(3),
38     &  shp2(7,8),vold(0:mi(2),*),xi,et,weight,voldprev(0:mi(2),*),
39     &  xlocal20(3,9,6),xlocal4(3,1,4),xlocal10(3,3,4),xlocal6(3,1,5),
40     &  xlocal15(3,4,5),xlocal8(3,4,6),xlocal8r(3,1,6),
41     &  disp(3),xload(2,*),delexternalwork
42!
43      include "gauss.f"
44      include "xlocal.f"
45!
46      data ifaceq /4,3,2,1,11,10,9,12,
47     &            5,6,7,8,13,14,15,16,
48     &            1,2,6,5,9,18,13,17,
49     &            2,3,7,6,10,19,14,18,
50     &            3,4,8,7,11,20,15,19,
51     &            4,1,5,8,12,17,16,20/
52      data ifacet /1,3,2,7,6,5,
53     &             1,2,4,5,9,8,
54     &             2,3,4,6,10,9,
55     &             1,4,3,8,10,7/
56      data ifacew /1,3,2,9,8,7,0,0,
57     &             4,5,6,10,11,12,0,0,
58     &             1,2,5,4,7,14,10,13,
59     &             2,3,6,5,8,15,11,14,
60     &             4,6,3,1,12,15,9,13/
61      data iflag /3/
62!
63!           determining the set with the external faces
64!
65      ipos=index(externalfaces,' ')
66      externalfaces(ipos:ipos)='T'
67      do iset=1,nset
68         if(set(iset)(1:ipos).eq.externalfaces(1:ipos)) exit
69      enddo
70      externalfaces(ipos:ipos)=' '
71!
72      delexternalwork=0.d0
73!
74      do jj=istartset(iset),iendset(iset)
75!
76         jface=ialset(jj)
77!
78         nelem=int(jface/10.d0)
79c** start change 20191219
80         indexe=ipkon(nelem)
81         if(indexe.lt.0) cycle
82c** end chenge 20191219
83         ig=jface-10*nelem
84!
85!        determining the distributed load on the face
86!
87         ipres=0
88         call nident2(nelemload,nelem,nload,id)
89         do
90            if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit
91            if(sideload(id)(1:1).ne.'P') then
92               id=id-1
93               cycle
94            endif
95            igpres=ichar(sideload(id)(2:2))-48
96            if(ig.ne.igpres) then
97               id=id-1
98               cycle
99            else
100               ipres=1
101               pres=xload(1,id)
102               exit
103            endif
104         enddo
105c         if(ipres.eq.0) then
106c            write(*,*) '*ERROR in calcexternalwork: external face',
107c     &           ig,' of element ',nelem,' has no external load'
108c            call exit(201)
109c         endif
110!
111c         indexe=ipkon(nelem)
112         lakonl=lakon(nelem)
113!
114!        determining the number of nodes in the face
115!
116         if(lakonl(4:4).eq.'2') then
117            nope=20
118            nopes=8
119         elseif(lakonl(4:4).eq.'8') then
120            nope=8
121            nopes=4
122         elseif(lakonl(4:5).eq.'10') then
123            nope=10
124            nopes=6
125         elseif(lakonl(4:4).eq.'4') then
126            nope=4
127            nopes=3
128         elseif(lakonl(4:5).eq.'15') then
129            nope=15
130         elseif(lakonl(4:4).eq.'6') then
131            nope=6
132         endif
133!
134!        determining the number of integration points in the face
135!
136         if(lakonl(4:5).eq.'8R') then
137            mint2d=1
138         elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
139     &           then
140            if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'S').or.
141     &           (lakonl(7:7).eq.'E')) then
142               mint2d=2
143            else
144               mint2d=4
145            endif
146         elseif(lakonl(4:4).eq.'2') then
147            mint2d=9
148         elseif(lakonl(4:5).eq.'10') then
149            mint2d=3
150         elseif(lakonl(4:4).eq.'4') then
151            mint2d=1
152         endif
153!
154!        local topology
155!
156         do i=1,nope
157            konl(i)=kon(indexe+i)
158         enddo
159!
160!        treatment of wedge faces
161!
162         if(lakonl(4:4).eq.'6') then
163            mint2d=1
164            if(ig.le.2) then
165               nopes=3
166            else
167               nopes=4
168            endif
169         endif
170         if(lakonl(4:5).eq.'15') then
171            if(ig.le.2) then
172               mint2d=3
173               nopes=6
174            else
175               mint2d=4
176               nopes=8
177            endif
178         endif
179
180!
181!        coordinates and differential displacements of the nodes belonging
182!        to the face face
183!
184         if((nope.eq.20).or.(nope.eq.8)) then
185            do i=1,nopes
186               do j=1,3
187                  voldl2(j,i)=vold(j,konl(ifaceq(i,ig)))
188                  xl2(j,i)=co(j,konl(ifaceq(i,ig)))
189     &                 +voldl2(j,i)
190                  voldl2(j,i)=voldl2(j,i)-voldprev(j,konl(ifaceq(i,ig)))
191               enddo
192            enddo
193         elseif((nope.eq.10).or.(nope.eq.4)) then
194            do i=1,nopes
195               do j=1,3
196                  voldl2(j,i)=vold(j,konl(ifacet(i,ig)))
197                  xl2(j,i)=co(j,konl(ifacet(i,ig)))
198     &                 +voldl2(j,i)
199                 voldl2(j,i)=voldl2(j,i)-voldprev(j,konl(ifacet(i,ig)))
200               enddo
201            enddo
202         else
203            do i=1,nopes
204               do j=1,3
205                  voldl2(j,i)=vold(j,konl(ifacew(i,ig)))
206                  xl2(j,i)=co(j,konl(ifacew(i,ig)))
207     &                 +voldl2(j,i)
208                  voldl2(j,i)=voldl2(j,i)-voldprev(j,konl(ifacew(i,ig)))
209               enddo
210            enddo
211         endif
212!
213         do i=1,mint2d
214!
215!          local coordinates of the surface integration
216!           point within the surface local coordinate system
217!
218            if((lakonl(4:5).eq.'8R').or.
219     &           ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
220               xi=gauss2d1(1,i)
221               et=gauss2d1(2,i)
222               weight=weight2d1(i)
223            elseif((lakonl(4:4).eq.'8').or.
224     &              (lakonl(4:6).eq.'20R').or.
225     &              ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
226               xi=gauss2d2(1,i)
227               et=gauss2d2(2,i)
228               weight=weight2d2(i)
229            elseif(lakonl(4:4).eq.'2') then
230               xi=gauss2d3(1,i)
231               et=gauss2d3(2,i)
232               weight=weight2d3(i)
233            elseif((lakonl(4:5).eq.'10').or.
234     &              ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
235               xi=gauss2d5(1,i)
236               et=gauss2d5(2,i)
237               weight=weight2d5(i)
238            elseif((lakonl(4:4).eq.'4').or.
239     &              ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
240               xi=gauss2d4(1,i)
241               et=gauss2d4(2,i)
242               weight=weight2d4(i)
243            endif
244!
245!           local surface normal
246!
247            if(nopes.eq.8) then
248               call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
249            elseif(nopes.eq.4) then
250               call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
251            elseif(nopes.eq.6) then
252               call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
253            else
254               call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
255            endif
256!
257!           differential displacements at the integration point
258!
259            do j1=1,3
260               disp(j1)=0.d0
261               do i1=1,nopes
262                  disp(j1)=disp(j1)+shp2(4,i1)*voldl2(j1,i1)
263               enddo
264            enddo
265!
266!           total external work (pressure is a negative load)
267!
268            delexternalwork=delexternalwork-pres*weight*
269     &        (disp(1)*xsj2(1)+disp(2)*xsj2(2)+disp(3)*xsj2(3))
270         enddo
271      enddo
272!
273      return
274      end
275