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 distributesens(istartdesi,ialdesi,ipkon,lakon,
20     &   ipoface,ndesi,nodedesi,nodface,kon,co,dgdx,nobject,
21     &   weightformgrad,nodedesiinv,iregion,objectset,
22     &   dgdxglob,nk,physcon,nobjectstart)
23!
24!     calculation of the relationship between facial distributions
25!     and corresponding nodal values
26!
27!     transforming the nodal sensitivities into facially
28!     distributed sensitivities
29!
30      implicit none
31!
32      character*81 objectset(5,*)
33      character*8 lakon(*)
34!
35      integer idesvar,j,k,kk,l,m,m1,n,istartdesi(*),
36     &   ielem,ipoface(*),nodface(5,*),ndesi,nodedesi(*),nope,
37     &   ialdesi(*),nopes,indexe,ifour,indexel,
38     &   mint3d,mint2d,ipkon(*),konl(26),iflag,ifaceq(8,6),
39     &   ifacet(6,4),ifacew1(4,5),ifacew2(8,5),kon(*),
40     &   nodes1(8),ifacel,iobject,nodes2(8),indexs,nk,
41     &   ithree,i,node,ifacq(2,3,20),ifact(2,3,10),ifacw(2,3,15),
42     &   nopedesi,nnodes,nodedesiinv(*),nobject,iregion,iactnode,
43     &   nobjectstart
44!
45      real*8 xi,et,weight,xl(3,9),xs(3,2),xsj(3),shp(7,9),
46     &   co(3,*),xsjj,weightformgrad(ndesi),dgdx(ndesi,nobject),
47     &   dgdxglob(2,nk,nobject),physcon(*)
48!
49!
50!
51!     flag for shape functions
52!
53      data iflag /3/
54      data indexel /0/
55!
56      save indexel
57!
58      include "gauss.f"
59!
60      data ifaceq /4,3,2,1,11,10,9,12,
61     &            5,6,7,8,13,14,15,16,
62     &            1,2,6,5,9,18,13,17,
63     &            2,3,7,6,10,19,14,18,
64     &            3,4,8,7,11,20,15,19,
65     &            4,1,5,8,12,17,16,20/
66!
67!     nodes per face for tet elements
68!
69      data ifacet /1,3,2,7,6,5,
70     &             1,2,4,5,9,8,
71     &             2,3,4,6,10,9,
72     &             1,4,3,8,10,7/
73!
74!     nodes per face for linear wedge elements
75!
76      data ifacew1 /1,3,2,0,
77     &             4,5,6,0,
78     &             1,2,5,4,
79     &             2,3,6,5,
80     &             3,1,4,6/
81!
82!     nodes per face for quadratic wedge elements
83!
84      data ifacew2 /1,3,2,9,8,7,0,0,
85     &             4,5,6,10,11,12,0,0,
86     &             1,2,5,4,7,14,10,13,
87     &             2,3,6,5,8,15,11,14,
88     &             3,1,4,6,9,13,12,15/
89!
90!     ifacq indicates for each local node number to which
91!     faces this node belongs and the internal number
92!     within this face. e.g., local node number 1 belongs to
93!     face 1 with internal number 4, face 3 with internal number
94!     1 and face 6 with internal number 2 (first line underneath)
95!
96!     ifacq, ifact and ifacw can be derived in a unique way from
97!     ifaceq, ifacet and ifacew2
98!
99      data ifacq /1,4,3,1,6,2,
100     &1,3,3,2,4,1,
101     &1,2,4,2,5,1,
102     &1,1,5,2,6,1,
103     &2,1,3,4,6,3,
104     &2,2,3,3,4,4,
105     &2,3,4,3,5,4,
106     &2,4,5,3,6,4,
107     &1,7,3,5,0,0,
108     &1,6,4,5,0,0,
109     &1,5,5,5,0,0,
110     &1,8,6,5,0,0,
111     &2,5,3,7,0,0,
112     &2,6,4,7,0,0,
113     &2,7,5,7,0,0,
114     &2,8,6,7,0,0,
115     &3,8,6,6,0,0,
116     &3,6,4,8,0,0,
117     &4,6,5,8,0,0,
118     &5,6,6,8,0,0/
119!
120      data ifact /1,1,2,1,4,1,
121     &1,3,2,2,3,1,
122     &1,2,3,2,4,3,
123     &2,3,3,3,4,2,
124     &1,6,2,4,0,0,
125     &1,5,3,4,0,0,
126     &1,4,4,6,0,0,
127     &2,6,4,4,0,0,
128     &2,5,3,6,0,0,
129     &3,5,4,5,0,0/
130!
131      data ifacw /1,1,3,1,5,2,
132     &1,3,3,2,4,1,
133     &1,2,4,2,5,1,
134     &2,1,3,4,5,3,
135     &2,2,3,3,4,4,
136     &2,3,4,3,5,4,
137     &1,6,3,5,0,0,
138     &1,5,4,5,0,0,
139     &1,4,5,5,0,0,
140     &2,4,3,7,0,0,
141     &2,5,4,7,0,0,
142     &2,6,5,7,0,0,
143     &3,8,5,6,0,0,
144     &3,6,4,8,0,0,
145     &4,6,5,8,0,0/
146!
147!     flag for shape functions
148!
149      ifour=4
150      ithree=3
151!
152!
153!     Loop over all designvariables
154!
155      do idesvar=1,ndesi
156!
157         node=nodedesi(idesvar)
158!
159!     Loop over all elements belonging to the designvariable
160!
161         do j=istartdesi(idesvar),istartdesi(idesvar+1)-1
162            ielem=ialdesi(j)
163
164            indexe=ipkon(ielem)
165!
166!     mint2d: # of integration points on the surface
167!     nopes: # of nodes in the surface
168!     nope: # of nodes in the element
169!
170            if(lakon(ielem)(4:5).eq.'8R') then
171               mint2d=1
172               nopes=4
173               nope=8
174               nopedesi=3
175            elseif(lakon(ielem)(4:4).eq.'8') then
176               mint2d=4
177               nopes=4
178               nope=8
179               nopedesi=3
180            elseif(lakon(ielem)(4:5).eq.'20') then
181               mint2d=9
182               nopes=8
183               nope=20
184               nopedesi=5
185            elseif(lakon(ielem)(4:5).eq.'10') then
186               mint2d=3
187               nopes=6
188               nope=10
189c               nopedesi=3
190               nopedesi=4
191            elseif(lakon(ielem)(4:4).eq.'4') then
192               mint2d=1
193               nopes=3
194               nope=4
195               nopedesi=3
196            elseif(lakon(ielem)(4:4).eq.'6') then
197               mint2d=1
198               nopes=3
199               nope=6
200               nopedesi=3
201            elseif(lakon(ielem)(4:5).eq.'15') then
202               mint2d=3
203               nopes=6
204               nope=15
205c               nopedesi=3
206            else
207               exit
208            endif
209            if(iregion.eq.0) nopedesi=0
210!
211!     actual position of the nodes belonging to the
212!     master surface
213!
214            do l=1,nope
215               konl(l)=kon(indexe+l)
216            enddo
217!
218            do i=1,nope
219               if(kon(indexe+i).eq.node) exit
220            enddo
221!
222!     Loop over all faces of the actual element
223!
224            loop1: do m=1,3
225!
226!     hexahedral elements
227!
228               if((nope.eq.20).or.(nope.eq.8)) then
229                  k=ifacq(1,m,i)
230                  m1=ifacq(2,m,i)
231                  do l=1,nopes
232                     nodes1(l)=konl(ifaceq(l,k))
233                     nodes2(l)=nodes1(l)
234                  enddo
235                  call insertsorti(nodes1,ifour)
236c                  call isortii(nodes1,iaux,ifour,kflag)
237!
238!     tetrahedral element
239!
240               elseif((nope.eq.10).or.(nope.eq.4)) then
241                  k=ifact(1,m,i)
242                  m1=ifact(2,m,i)
243                  do l=1,nopes
244                     nodes1(l)=konl(ifacet(l,k))
245                     nodes2(l)=nodes1(l)
246                  enddo
247                  call insertsorti(nodes1,ithree)
248c                  call isortii(nodes1,iaux,ithree,kflag)
249                  nodes1(4)=0
250!
251!     wedge element
252!
253               elseif(nope.eq.15) then
254                  k=ifacw(1,m,i)
255                  m1=ifacw(2,m,i)
256                  if(k.le.2) then
257                     nopes=6
258                     nopedesi=4
259                     mint3d=3
260                     do l=1,nopes
261                        nodes1(l)=konl(ifacew2(l,k))
262                        nodes2(l)=nodes1(l)
263                     enddo
264                     call insertsorti(nodes1,ithree)
265c                     call isortii(nodes1,iaux,ithree,kflag)
266                     nodes1(4)=0
267                  else
268                     nopes=8
269                     nopedesi=5
270                     mint3d=4
271                     do l=1,nopes
272                        nodes1(l)=konl(ifacew2(l,k))
273                        nodes2(l)=nodes1(l)
274                     enddo
275                     call insertsorti(nodes1,ithree)
276c                     call isortii(nodes1,iaux,ithree,kflag)
277                  endif
278               else
279                  k=ifacw(1,m,i)
280                  m1=ifacw(2,m,i)
281                  if(k.le.2) then
282                     nopes=3
283                     mint3d=1
284                     do l=1,nopes
285                        nodes1(l)=konl(ifacew1(l,k))
286                        nodes2(l)=nodes1(l)
287                     enddo
288                     call insertsorti(nodes1,ithree)
289c                     call isortii(nodes1,iaux,ithree,kflag)
290                     nodes1(4)=0
291                  else
292                     nopes=4
293                     mint3d=2
294                     do l=1,nopes
295                        nodes1(l)=konl(ifacew1(l,k))
296                        nodes2(l)=nodes1(l)
297                     enddo
298                     call insertsorti(nodes1,ithree)
299c                     call isortii(nodes1,iaux,ithree,kflag)
300                  endif
301               endif
302               if(k.eq.0) exit
303!
304!              Find external element face
305!
306               indexs=ipoface(nodes1(1))
307               do
308                  if(indexs.eq.0) cycle loop1
309                  if((nodface(1,indexs).eq.nodes1(2)).and.
310     &                 (nodface(2,indexs).eq.nodes1(3))) then
311!
312!                    Check if sufficient designnodes on that surface
313!
314                     nnodes=0
315                     do n=1,nopes
316                        if(nodedesiinv(nodes1(n)).eq.1) then
317                           nnodes=nnodes+1
318                        endif
319                     enddo
320!
321!                    Assign coordinates to nodes of surface
322!
323c                     write(*,*) 'formgradient',nnodes,nopedesi
324                     if(nnodes.ge.nopedesi) then
325                        if((nope.eq.20).or.(nope.eq.8)) then
326                           do l=1,nopes
327                              do n=1,3
328                                 ifacel=ifaceq(l,nodface(4,indexs))
329                                 xl(n,l)=co(n,konl(ifacel))
330                              enddo
331                           enddo
332                        elseif((nope.eq.10).or.(nope.eq.4)) then
333                           do l=1,nopes
334                              do n=1,3
335                                 ifacel=ifacet(l,nodface(4,indexs))
336                                 xl(n,l)=co(n,konl(ifacel))
337                              enddo
338                           enddo
339                        elseif(nope.eq.15) then
340                           do l=1,nopes
341                              do n=1,3
342                                 ifacel=ifacew2(l,nodface(4,indexs))
343                                 xl(n,l)=co(n,konl(ifacel))
344                              enddo
345                           enddo
346                        elseif(nope.eq.6) then
347                           do l=1,nopes
348                              do n=1,3
349                                 ifacel=ifacew1(l,nodface(4,indexs))
350                                 xl(n,l)=co(n,konl(ifacel))
351                              enddo
352                           enddo
353                        endif
354!
355!                       Determine weighting factors
356!
357                        do kk=1,mint2d
358                           if(lakon(ielem)(4:5).eq.'20') then
359                              xi=gauss2d3(1,kk)
360                              et=gauss2d3(2,kk)
361                              weight=weight2d3(kk)
362                              call shape8q(xi,et,xl,xsj,
363     &                             xs,shp,iflag)
364                           elseif(lakon(ielem)(4:5).eq.'10') then
365                              xi=gauss2d5(1,kk)
366                              et=gauss2d5(2,kk)
367                              weight=weight2d5(kk)
368                              call shape6tri(xi,et,xl,xsj,
369     &                             xs,shp,iflag)
370                           elseif(lakon(ielem)(4:4).eq.'4') then
371                              xi=gauss2d4(1,kk)
372                              et=gauss2d4(2,kk)
373                              weight=weight2d4(kk)
374                              call shape3tri(xi,et,xl,xsj,
375     &                             xs,shp,iflag)
376                           elseif(lakon(ielem)(4:5).eq.'15') then
377                              if(k.le.2) then
378                                 xi=gauss2d5(1,kk)
379                                 et=gauss2d5(2,kk)
380                                 weight=weight2d5(kk)
381                                 call shape6tri(xi,et,xl,xsj,
382     &                                xs,shp,iflag)
383                              else
384                                 xi=gauss2d3(1,kk)
385                                 et=gauss2d3(2,kk)
386                                 weight=weight2d3(kk)
387                                 call shape8q(xi,et,xl,xsj,
388     &                                xs,shp,iflag)
389                              endif
390                           elseif(lakon(ielem)(4:5).eq.'8R') then
391                              xi=gauss2d1(1,kk)
392                              et=gauss2d1(2,kk)
393                              weight=weight2d1(kk)
394                              call shape4q(xi,et,xl,xsj,
395     &                             xs,shp,iflag)
396                           elseif(lakon(ielem)(4:4).eq.'8') then
397                              xi=gauss2d2(1,kk)
398                              et=gauss2d2(2,kk)
399                              weight=weight2d2(kk)
400                              call shape4q(xi,et,xl,xsj,
401     &                             xs,shp,iflag)
402                           elseif(lakon(ielem)(4:4).eq.'6') then
403                              if(k.le.2) then
404                                 xi=gauss2d4(1,kk)
405                                 et=gauss2d4(2,kk)
406                                 weight=weight2d4(kk)
407                                 call shape3tri(xi,et,xl,xsj,
408     &                                xs,shp,iflag)
409                               else
410                                 xi=gauss2d2(1,kk)
411                                 et=gauss2d2(2,kk)
412                                 weight=weight2d2(kk)
413                                 call shape4q(xi,et,xl,xsj,
414     &                                xs,shp,iflag)
415                              endif
416                           endif
417!
418!                          Calculate Jacobian determinant
419!
420                           xsjj=dsqrt(xsj(1)**2+xsj(2)**2+
421     &                          xsj(3)**2)
422!
423!                          Evaluate Shape functions for weightmatrix
424!
425                           weightformgrad(idesvar)=weightformgrad
426     &                          (idesvar)+weight*shp(4,m1)*xsjj
427!                           write(*,*) 'formgradient ',idesvar,j,m,kk,
428!     &                             weightformgrad(idesvar)
429                        enddo
430                     endif
431                     exit
432                  endif
433                  indexs=nodface(5,indexs)
434               enddo            ! find external element faces
435!
436            enddo loop1         ! Loop over all faces of the actual element
437!
438         enddo                  ! Loop over all elements belonging to the designvariable
439!
440      enddo                     ! Loop over all designvariables
441!
442!     Change sign of sensitivity of objective function in case of
443!     a maximization problem
444!
445      if(physcon(11).le.0.d0) then
446         if(objectset(2,1)(17:19).eq.'MAX') then
447            do idesvar=1,ndesi
448               dgdx(idesvar,1)=-dgdx(idesvar,1)
449            enddo
450         endif
451      endif
452!
453!     Scaling of sensitivities
454!
455      do idesvar=1,ndesi
456         do iobject=1+nobjectstart,nobject
457            if(objectset(1,iobject)(4:13).eq.'MEMBERSIZE') cycle
458            if(objectset(1,iobject)(1:9).eq.'FIXGROWTH') cycle
459            if(objectset(1,iobject)(1:12).eq.'FIXSHRINKAGE') cycle
460            if(weightformgrad(idesvar).gt.0.d0) then
461               dgdx(idesvar,iobject)=dgdx(idesvar,iobject)
462     &             /weightformgrad(idesvar)
463               iactnode=nodedesi(idesvar)
464               dgdxglob(1,iactnode,iobject)=dgdx(idesvar,iobject)
465            else
466               iactnode=nodedesi(idesvar)
467               dgdxglob(1,iactnode,iobject)=dgdx(idesvar,iobject)
468            endif
469         enddo
470      enddo
471!
472      return
473      end
474
475