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 getdesiinfo3d(set,istartset,iendset,ialset,nset,
20     &     mi,nactdof,ndesi,nodedesi,ntie,tieset,itmp,nmpc,nodempc,
21     &     ipompc,nodedesiinv,iponoel,inoel,lakon,ipkon,kon,iregion,
22     &     ipoface,nodface,nk)
23!
24!     storing the design variables in nodedesi
25!     marking which nodes are design variables in nodedesiinv
26!
27!     a node is a design variable if:
28!     1) it belongs to the design variable set AND
29!     2) not all dofs in the node are defined by SPC's AND
30!     3) no MPC is applied to any of its dofs AND
31!     4) it belongs to at least one face whose number of
32!     design variables exceeds half its nodes
33!
34      implicit none
35!
36      character*8 lakon(*)
37!
38      character*81 setname
39      character*81 set(*)
40      character*81 tieset(3,*)
41!
42      integer mi(*),istartset(*),iendset(*),ialset(*),ndesi,
43     &     node,nodedesi(*),nset,ntie,i,j,k,l,m,nmpc,nodempc(3,*),
44     &     nactdof(0:mi(2),*),itmp(*),ntmp,index,id,ipompc(*),
45     &     nodedesiinv(*),iponoel(*),inoel(2,*),nelem,nope,nopedesi,
46     &     ipkon(*),nnodes,kon(*),iregion,konl(26),iaux,kflag,
47     &     ipoface(*),nodface(5,*),jfacem,nopesurf(9),ifaceq(8,6),
48     &     ifacet(6,4),ifacew1(4,5),ifacew2(8,5),nopem,nk,iwrite
49!
50      setname(1:1)=' '
51      ndesi=0
52!
53!     nodes per face for hex elements
54!
55      data ifaceq /4,3,2,1,11,10,9,12,
56     &     5,6,7,8,13,14,15,16,
57     &     1,2,6,5,9,18,13,17,
58     &     2,3,7,6,10,19,14,18,
59     &     3,4,8,7,11,20,15,19,
60     &     4,1,5,8,12,17,16,20/
61!
62!     nodes per face for tet elements
63!
64      data ifacet /1,3,2,7,6,5,
65     &     1,2,4,5,9,8,
66     &     2,3,4,6,10,9,
67     &     1,4,3,8,10,7/
68!
69!     nodes per face for linear wedge elements
70!
71      data ifacew1 /1,3,2,0,
72     &     4,5,6,0,
73     &     1,2,5,4,
74     &     2,3,6,5,
75     &     3,1,4,6/
76!
77!     nodes per face for quadratic wedge elements
78!
79      data ifacew2 /1,3,2,9,8,7,0,0,
80     &     4,5,6,10,11,12,0,0,
81     &     1,2,5,4,7,14,10,13,
82     &     2,3,6,5,8,15,11,14,
83     &     3,1,4,6,9,13,12,15/
84!
85!     Search for the set name of the set with the design variables
86!
87      do i=1,ntie
88        if(tieset(1,i)(81:81).eq.'D') then
89          setname=tieset(2,i)
90        endif
91      enddo
92!
93!     Check for the existence of the name
94!
95      if(setname(1:1).eq.' ') then
96        write(*,*) '*ERROR in getdesiinfo: name of node set '
97        write(*,*) '  has not yet been defined. '
98        call exit(201)
99      endif
100!
101!     catalogue all nodes (dependent and independent) which
102!     belong to MPC's and sort them in increasing order
103!
104      ntmp=0
105      do i=1,nmpc
106        index=ipompc(i)
107        do
108          if(index.eq.0) exit
109          node=nodempc(1,index)
110          call nident(itmp,node,ntmp,id)
111          if(id.gt.0) then
112            if(itmp(id).eq.node) then
113              index=nodempc(3,index)
114              cycle
115            endif
116          endif
117          ntmp=ntmp+1
118          do j=ntmp,id+2,-1
119            itmp(j)=itmp(j-1)
120          enddo
121          itmp(id+1)=node
122          index=nodempc(3,index)
123        enddo
124      enddo
125!
126!     opening a file to store the nodes which are rejected as
127!     design variables
128!
129      iwrite=0
130      open(40,file='WarnNodeDesignReject.nam',status='unknown')
131      write(40,*) '*NSET,NSET=WarnNodeDesignReject'
132!
133!     Search the name of the node set in "set(i)" and
134!     assign the nodes of the set to the appropriate variables
135!
136c      do i=1,nset
137c        if(setname.eq.set(i)) then
138      call cident81(set,setname,nset,i)
139      if(i.gt.0) then
140        if(setname.eq.set(i)) then
141          loop1: do j=istartset(i),iendset(i)
142            if(ialset(j).gt.0) then
143              node=ialset(j)
144!
145!     check for SPC-constraints: if a node is constrained in
146!     all dofs it is removed from the design node set
147!
148              do l=1,3
149                if(nactdof(l,node).gt.0) exit
150                if(l.eq.3) then
151                  write(*,*) '*WARNING in getdesiinfo:'
152                  write(*,*) '         node ',node,' has no'
153                  write(*,*) '         active dofs and'
154                  write(*,*) '         is removed from the set'
155                  write(*,*) '         of design variables'
156                  write(40,*) node
157                  iwrite=1
158                  cycle loop1
159                endif
160              enddo
161!
162!     check for MPC-constraints
163!
164              call nident(itmp,node,ntmp,id)
165              if(id.gt.0) then
166                if(itmp(id).eq.node) then
167                  write(*,*) '*WARNING in getdesiinfo:'
168                  write(*,*) '       node ',node,' is subject'
169                  write(*,*) '       to MPC-constraints and'
170                  write(*,*) '       is removed from the set'
171                  write(*,*) '       of design variables'
172                  write(40,*) node
173                  iwrite=1
174                  cycle loop1
175                endif
176              endif
177!
178              ndesi=ndesi+1
179              nodedesi(ndesi)=node
180            else
181              k=ialset(j-2)
182              loop2: do
183                k=k-ialset(j)
184                if(k.ge.ialset(j-1)) exit
185!
186!     check for SPC-constraints: if a node is constrained in
187!     all dofs it is removed from the design node set
188!
189                do l=1,3
190                  if(nactdof(l,k).gt.0) exit
191                  if(l.eq.3) then
192                    write(*,*) '*WARNING in getdesiinfo:'
193                    write(*,*) '         node ',k,' has no'
194                    write(*,*) '         active dofs and'
195                    write(*,*) '         is removed from the set'
196                    write(*,*) '         of design variables'
197                    write(40,*) k
198                    iwrite=1
199                    cycle loop2
200                  endif
201                enddo
202!
203!     check for MPC-constraints
204!
205                call nident(itmp,k,ntmp,id)
206                if(id.gt.0) then
207                  if(itmp(id).eq.k) then
208                    write(*,*) '*WARNING in getdesiinfo:'
209                    write(*,*) '   node ',k,' is subject'
210                    write(*,*) '   to MPC-constraints and'
211                    write(*,*) '   is removed from the set'
212                    write(*,*) '   of design variables'
213                    write(40,*) k
214                    iwrite=1
215                    cycle loop2
216                  endif
217                endif
218!
219                ndesi=ndesi+1
220                nodedesi(ndesi)=k
221              enddo loop2
222            endif
223          enddo loop1
224        endif
225      endif
226c        endif
227c      enddo
228!
229!     creating field nodedesiinv indicating for each node whether
230!     it is a design variable or not
231!
232      do i=1,ndesi
233        index=nodedesi(i)
234        nodedesiinv(index)=-1
235      enddo
236!
237      kflag=1
238      call isortii(nodedesi,iaux,ndesi,kflag)
239!
240!     A design node is also removed from nodedesi if it does not
241!     belong to a face whose number of design variables exceeds half
242!     of its nodes
243!
244      do i=1,nk
245        node=i
246        if(ipoface(node).eq.0) cycle
247        index=ipoface(node)
248        do
249          nelem=nodface(3,index)
250          jfacem=nodface(4,index)
251!
252          if(lakon(nelem)(4:4).eq.'8') then
253            nope=8
254            nopedesi=3
255            nopem=4
256          elseif(lakon(nelem)(4:5).eq.'20') then
257            nope=20
258            nopedesi=5
259            nopem=8
260          elseif(lakon(nelem)(4:5).eq.'10') then
261            nope=10
262            nopedesi=4
263            nopem=6
264          elseif(lakon(nelem)(4:4).eq.'4') then
265            nope=4
266            nopedesi=3
267            nopem=3
268          elseif(lakon(nelem)(4:4).eq.'6') then
269            nope=6
270            if(jfacem.le.2) then
271              nopem=3
272              nopedesi=3
273            else
274              nopem=4
275              nopedesi=3
276            endif
277          elseif(lakon(nelem)(4:5).eq.'15') then
278            nope=15
279            if(jfacem.le.2) then
280              nopem=6
281              nopedesi=4
282            else
283              nopem=8
284              nopedesi=5
285            endif
286          endif
287          if(iregion.eq.0) nopedesi=0
288!
289!     actual position of the nodes belonging to the
290!     master surface
291!
292          do k=1,nope
293            konl(k)=kon(ipkon(nelem)+k)
294          enddo
295!
296          if((nope.eq.20).or.(nope.eq.8)) then
297            do m=1,nopem
298              nopesurf(m)=konl(ifaceq(m,jfacem))
299            enddo
300          elseif((nope.eq.10).or.(nope.eq.4)) then
301            do m=1,nopem
302              nopesurf(m)=konl(ifacet(m,jfacem))
303            enddo
304          elseif(nope.eq.15) then
305            do m=1,nopem
306              nopesurf(m)=konl(ifacew2(m,jfacem))
307            enddo
308          else
309            do m=1,nopem
310              nopesurf(m)=konl(ifacew1(m,jfacem))
311            enddo
312          endif
313!
314!     sum up how many designvariables are on that surface
315!
316          nnodes=0
317          do m=1,nopem
318            if(nodedesiinv(nopesurf(m)).ne.0) then
319              nnodes=nnodes+1
320            endif
321          enddo
322!
323          if(nnodes.ge.nopedesi) then
324            do m=1,nopem
325              if(nodedesiinv(nopesurf(m)).eq.-1) then
326                nodedesiinv(nopesurf(m))=1
327              endif
328            enddo
329          endif
330          index=nodface(5,index)
331          if(index.eq.0) exit
332        enddo
333      enddo
334!
335!     if node i in nodedesi(i) is -1 --> delete node i from
336!     set of designvariables
337!
338      do i=1,nk
339        if(nodedesiinv(i).eq.-1) then
340!
341          write(*,*) '*WARNING in getdesiinfo:'
342          write(*,*) '          node ',i,' is removed'
343          write(*,*) '          from the set of design'
344          write(*,*) '          variables as not sufficient '
345          write(*,*) '          other variables are on the  '
346          write(*,*) '          surrounding element faces  '
347          write(40,*) i
348          iwrite=1
349!
350          nodedesiinv(i)=0
351          call nident(nodedesi,i,ndesi,id)
352          do k=id+1,ndesi
353            nodedesi(k-1)=nodedesi(k)
354          enddo
355          ndesi=ndesi-1
356        endif
357      enddo
358!
359      if(iwrite.eq.1) then
360        write(*,*) '*INFO in getdesiinfo:'
361        write(*,*) '      rejected design nodes are stored in'
362        write(*,*) '      file WarnNodeDesignReject.nam'
363        write(*,*) '      This file can be loaded into'
364        write(*,*) '      an active cgx-session by typing'
365        write(*,*)
366     &       '      read WarnNodeDesignReject.nam inp'
367        write(*,*)
368        close(40)
369      else
370        close(40,status='delete')
371      endif
372!
373      return
374      end
375