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 patch(iterms,node,sti,scpav,mi,kon,ipkon,
20     & ipoints,members,linpatch,co,lakon,iavflag)
21!
22!     computes the smoothed nodal stresses for an element patch
23!
24!     author: Sascha Merz
25!
26      implicit none
27!
28      integer i,j,nope,mint3d,indexe,ielem,kon(*),ipkon(*),
29     & iterms,k,iflag,nrhs,info,node,ipnt,
30     & mi(*),irow,ipoints,members(*),linpatch,ielidx,iavflag
31!
32      real*8 xl(3,20),co(3,*),shp(4,20),
33     & pgauss(3),pol(20),pntdist,w,scpav(6,*),xsj,
34     & sti(6,mi(1),*),xi,et,ze,rv1(ipoints),pp(ipoints,iterms),
35     & pdat(ipoints,6),pfit(6),pwrk(iterms),pre(ipoints,iterms),
36     & z(ipoints,ipoints)
37!
38      real*8 tmpstr(6),gauss3d5e(3,4)
39!
40      character*8 lakon(*)
41!
42      logical matu,matv
43!
44      include 'gauss.f'
45!
46!     initialize
47!
48!     iavflag: if 1, then build average of integration point values
49!
50      if(iavflag.eq.0) then
51         do i=1,6
52            pfit(i)=0.d0
53         enddo
54         iflag=1
55!
56!        irow: row number of the rectangular matrix of the overdetermined
57!        system of equations
58!
59         irow=0
60!
61!        loop over patch elements
62!        linpatch: number of elements in patch
63!
64         do ielidx=1,linpatch
65!
66            ielem=members(ielidx)
67            if(lakon(ielem)(1:5).eq.'C3D20') then
68!
69!              nope: nodes per element
70!
71               nope=20
72               if(lakon(ielem)(6:6).eq.'R') then
73                  mint3d=8
74               else
75                  mint3d=27
76               endif
77            elseif(lakon(ielem)(1:4).eq.'C3D8') then
78               nope=8
79               if(lakon(ielem)(5:5).eq.'R') then
80                  mint3d=1
81               else
82                  mint3d=8
83               endif
84            elseif(lakon(ielem)(1:5).eq.'C3D10') then
85               nope=10
86               mint3d=4
87            endif
88!
89            indexe=ipkon(ielem)
90!
91!           coordinates of the nodes belonging to the element
92!
93            do j=1,nope
94               do k=1,3
95                  xl(k,j)=co(k,kon(indexe+j))
96               enddo
97            enddo
98!
99!           loop over the integration points in one element
100!
101            do ipnt=1,mint3d
102               irow=irow+1
103               if(nope.eq.10) then
104                  xi=gauss3d5(1,ipnt)
105                  et=gauss3d5(2,ipnt)
106                  ze=gauss3d5(3,ipnt)
107                  call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
108               elseif(nope.eq.20) then
109                  if(mint3d.eq.8) then
110                     xi=gauss3d2(1,ipnt)
111                     et=gauss3d2(2,ipnt)
112                     ze=gauss3d2(3,ipnt)
113                  elseif(mint3d.eq.27) then
114                     xi=gauss3d3(1,ipnt)
115                     et=gauss3d3(2,ipnt)
116                     ze=gauss3d3(3,ipnt)
117                  endif
118                  call shape20h(xi,et,ze,xl,xsj,shp,iflag)
119               elseif(nope.eq.8) then
120                  if(mint3d.eq.1) then
121                     xi=gauss3d1(1,ipnt)
122                     et=gauss3d1(2,ipnt)
123                     ze=gauss3d1(3,ipnt)
124                  elseif(mint3d.eq.8) then
125                     xi=gauss3d2(1,ipnt)
126                     et=gauss3d2(2,ipnt)
127                     ze=gauss3d2(3,ipnt)
128                  endif
129                  call shape8h(xi,et,ze,xl,xsj,shp,iflag)
130               endif
131!
132!              in pgauss the global coordinates of the integration
133!              point are saved.
134!
135               do j=1,3
136                  pgauss(j)=0.d0
137                  do k=1,nope
138                     pgauss(j)=pgauss(j)+shp(4,k)*xl(j,k)
139                  enddo
140!
141!                 the origin of the coordinate system is moved to the
142!                 evaluated node for higher numerical stability
143!
144                  pgauss(j)=pgauss(j)-co(j,node)
145               enddo
146!
147!              evaluate patch polynomial for the integration point
148!
149               pol(1)=1.d0
150               pol(2)=pgauss(1)
151               pol(3)=pgauss(2)
152               pol(4)=pgauss(3)
153               pol(5)=pgauss(1)*pgauss(2)
154               pol(6)=pgauss(1)*pgauss(3)
155               pol(7)=pgauss(2)*pgauss(3)
156               pol(8)=pgauss(1)*pgauss(1)
157               pol(9)=pgauss(2)*pgauss(2)
158               pol(10)=pgauss(3)*pgauss(3)
159               pol(11)=pgauss(1)*pgauss(2)*pgauss(3)
160               if(iterms.gt.11) then
161                  pol(12)=pgauss(1)*pgauss(1)*pgauss(2)
162                  pol(13)=pgauss(1)*pgauss(2)*pgauss(2)
163                  pol(14)=pgauss(1)*pgauss(1)*pgauss(3)
164                  pol(15)=pgauss(1)*pgauss(3)*pgauss(3)
165                  pol(16)=pgauss(2)*pgauss(2)*pgauss(3)
166                  pol(17)=pgauss(2)*pgauss(3)*pgauss(3)
167                  if(iterms.gt.17) then
168                     pol(18)=pgauss(1)*pgauss(1)*pgauss(1)
169                     pol(19)=pgauss(2)*pgauss(2)*pgauss(2)
170                     pol(20)=pgauss(3)*pgauss(3)*pgauss(3)
171                  endif
172               endif
173!
174!              weighting for integration point
175!
176               pntdist=dsqrt(pgauss(1)*pgauss(1)
177     &              +pgauss(2)*pgauss(2)
178     &              +pgauss(3)*pgauss(3))
179               w=pntdist**(-1.5d0)
180!
181               do j=1,6
182                  pdat(irow,j)=sti(j,ipnt,ielem)*w
183               enddo
184!
185               do j=1,iterms
186                  pp(irow,j)=pol(j)*w
187               enddo
188            enddo
189         enddo
190!
191!        using singular value decomposition for the least squares fit
192!
193         matu=.false.
194         matv=.true.
195         nrhs=6
196!
197         call hybsvd(ipoints,ipoints,ipoints,ipoints,ipoints,
198     &        ipoints,iterms,pp,pwrk,matu,pp,matv,
199     &        pre,z,pdat,nrhs,info,rv1)
200         if(info.ne.0) then
201            write(*,*) '*ERROR in patch: Bad conditioned matrix,',
202     &           ' using average of sampling point values.'
203            iavflag=1
204         endif
205      endif
206!
207!     matrix multiplication. only the first value of the
208!     solution vector is needed. the singular values are manipulated
209!     to increase the numerical stability
210!
211      if(iavflag.eq.0) then
212         do j=1,iterms
213            if(pwrk(j).lt.1.d-22) then
214               pwrk(j)=0.d0
215            else
216               pwrk(j)=1.d0/pwrk(j)
217            endif
218         enddo
219         do j=1,iterms
220            pre(1,j)=pre(1,j)*pwrk(j)
221         enddo
222         do j=1,nrhs
223            do k=1,iterms
224               pfit(j)=pfit(j)+pre(1,k)*pdat(k,j)
225            enddo
226         enddo
227!
228!        pfit is an array containing the coefficients for the polynom
229!        for the six stress components
230!
231!        solution in the node
232!
233         do j=1,6
234            scpav(j,node)=scpav(j,node)+pfit(j)
235         enddo
236      endif
237!
238!     if there are not enough elements to fit a polynomial,
239!     build average value of the sampling point values
240!
241      if(iavflag.eq.1) then
242         do j=1,6
243            tmpstr(j)=0.d0
244         enddo
245         irow=0
246         do ielidx=1,linpatch
247            ielem=members(ielidx)
248            if(lakon(ielem)(1:5).eq.'C3D20') then
249               if(lakon(ielem)(6:6).eq.'R') then
250                  mint3d=8
251               else
252                  mint3d=27
253               endif
254            elseif(lakon(ielem)(1:5).eq.'C3D10') then
255               mint3d=4
256            elseif(lakon(ielem)(1:4).eq.'C3D8') then
257               if(lakon(ielem)(5:5).eq.'R') then
258                  mint3d=1
259               else
260                  mint3d=8
261               endif
262            endif
263            do ipnt=1,mint3d
264               irow=irow+1
265               do j=1,6
266                  tmpstr(j)=tmpstr(j)+sti(j,ipnt,ielem)
267               enddo
268            enddo
269         enddo
270         do j=1,6
271            scpav(j,node)=scpav(j,node)+tmpstr(j)/irow
272         enddo
273      endif
274!
275      return
276      end
277