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_robust(set,istartset,iendset,ialset,
20     &  nset,mi,nactdof,ndesi,nodedesi,ntie,tieset,itmp,nmpc,nodempc,
21     &  ipompc,nodedesiinv,iponoel,inoel,lakon,ipkon,kon,iregion,
22     &  ipoface,nodface,nk,irandomtype)
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 for the random field 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,ishift,
49     &  indexe,jface,konl2d(26),expandhex(20),expandwed(15),
50     &  idelta,idesi,irandomtype(*),iwrite
51!
52      setname(1:1)=' '
53      ndesi=0
54!
55!     nodes per face for hex elements
56!
57      data ifaceq /4,3,2,1,11,10,9,12,
58     &            5,6,7,8,13,14,15,16,
59     &            1,2,6,5,9,18,13,17,
60     &            2,3,7,6,10,19,14,18,
61     &            3,4,8,7,11,20,15,19,
62     &            4,1,5,8,12,17,16,20/
63!
64!     nodes per face for tet elements
65!
66      data ifacet /1,3,2,7,6,5,
67     &             1,2,4,5,9,8,
68     &             2,3,4,6,10,9,
69     &             1,4,3,8,10,7/
70!
71!     nodes per face for linear wedge elements
72!
73      data ifacew1 /1,3,2,0,
74     &             4,5,6,0,
75     &             1,2,5,4,
76     &             2,3,6,5,
77     &             3,1,4,6/
78!
79!     nodes per face for quadratic wedge elements
80!
81      data ifacew2 /1,3,2,9,8,7,0,0,
82     &             4,5,6,10,11,12,0,0,
83     &             1,2,5,4,7,14,10,13,
84     &             2,3,6,5,8,15,11,14,
85     &             3,1,4,6,9,13,12,15/
86!
87!     opening a file to store the nodes which are rejected as
88!     design variables
89!
90      iwrite=0
91      open(40,file='WarnNodeDesignReject.nam',status='unknown')
92      write(40,*) '*NSET,NSET=WarnNodeDesignReject'
93!
94!------------------------------------------------------------------------
95!     assign the nodes of the set to the appropriate variables
96!------------------------------------------------------------------------
97!
98      do i=1,nk
99         if(irandomtype(i).gt.0) then
100            ndesi=ndesi+1
101            nodedesi(ndesi)=i
102            nodedesiinv(i)=-2
103         endif
104      enddo
105!
106!------------------------------------------------------------------------
107!     Remove interior nodes from the set of design variables
108!------------------------------------------------------------------------
109!
110      do j=1,nk
111!
112         if(ipoface(j).eq.0) cycle
113         indexe=ipoface(j)
114!
115         do
116!
117            nelem=nodface(3,indexe)
118            jface=nodface(4,indexe)
119!
120            if((lakon(nelem)(7:7).eq.'A').or.
121     &         (lakon(nelem)(7:7).eq.'S').or.
122     &         (lakon(nelem)(7:7).eq.'E')) then
123!
124!              for plane stress/strain/axi only faces
125!              3 and higher are taken into account for the normal
126!
127               if(jface.le.2) then
128                  indexe=nodface(5,indexe)
129                  if(indexe.eq.0) then
130                     exit
131                  else
132                     cycle
133                  endif
134               endif
135            elseif(lakon(nelem)(7:7).eq.'L') then
136!
137!              for shells only faces 2 and lower
138!              taken into account for the normal
139!
140               if(jface.gt.2) then
141                  indexe=nodface(5,indexe)
142                  if(indexe.eq.0) then
143                     exit
144                  else
145                     cycle
146                  endif
147               endif
148            endif
149!
150!           nopem: # of nodes in the surface
151!           nope: # of nodes in the element
152!
153            if(lakon(nelem)(4:4).eq.'8') then
154               nopem=4
155               nope=8
156               nopedesi=3
157               ishift=8
158            elseif(lakon(nelem)(4:5).eq.'20') then
159               nopem=8
160               nope=20
161               nopedesi=5
162               ishift=20
163            elseif(lakon(nelem)(4:5).eq.'10') then
164               nopem=6
165               nope=10
166               nopedesi=4
167            elseif(lakon(nelem)(4:4).eq.'4') then
168               nopem=3
169               nope=4
170               nopedesi=3
171!
172!     treatment of wedge faces
173!
174            elseif(lakon(nelem)(4:4).eq.'6') then
175               nope=6
176               if(jface.le.2) then
177                  nopem=3
178               else
179                  nopem=4
180               endif
181               nopedesi=3
182               ishift=6
183            elseif(lakon(nelem)(4:5).eq.'15') then
184               nope=15
185               if(jface.le.2) then
186                  nopem=6
187                  nopedesi=4
188               else
189                  nopem=8
190                  nopedesi=5
191               endif
192               ishift=15
193            endif
194            if(iregion.eq.0) then
195         nopedesi=0
196      endif
197!
198!     actual position of the nodes belonging to the surface
199!
200            if((lakon(nelem)(7:7).eq.'A').or.
201     &         (lakon(nelem)(7:7).eq.'S').or.
202     &         (lakon(nelem)(7:7).eq.'E')) then
203               if((lakon(nelem)(4:5).eq.'20').or.
204     &            (lakon(nelem)(4:5).eq.'8 ')) then
205                  do k=1,nope
206                     konl(k)=kon(ipkon(nelem)+k)
207                     konl2d(k)=kon(ipkon(nelem)+ishift+expandhex(k))
208                  enddo
209               elseif((lakon(nelem)(4:5).eq.'15').or.
210     &                 (lakon(nelem)(4:5).eq.'6 ')) then
211                  do k=1,nope
212                     konl(k)=kon(ipkon(nelem)+k)
213                     konl2d(k)=kon(ipkon(nelem)+ishift+expandwed(k))
214                  enddo
215               endif
216            else
217               do k=1,nope
218                  konl(k)=kon(ipkon(nelem)+k)
219               enddo
220            endif
221!
222            if((nope.eq.20).or.(nope.eq.8)) then
223               do m=1,nopem
224                  nopesurf(m)=konl(ifaceq(m,jface))
225               enddo
226            elseif((nope.eq.10).or.(nope.eq.4)) then
227               do m=1,nopem
228                  nopesurf(m)=konl(ifacet(m,jface))
229               enddo
230            elseif(nope.eq.15) then
231               do m=1,nopem
232                  nopesurf(m)=konl(ifacew2(m,jface))
233               enddo
234            else
235               do m=1,nopem
236                  nopesurf(m)=konl(ifacew1(m,jface))
237               enddo
238            endif
239!
240!     sum up how many designvariables are on that surface
241!
242            nnodes=0
243            do m=1,nopem
244               if(nodedesiinv(nopesurf(m)).ne.0) then
245                  nnodes=nnodes+1
246               endif
247            enddo
248!
249!     set design variables located on external surfaces
250!     equal to 1 in nodedesiinv
251!     if sufficient design variables are defined on this surface
252!
253            nnodes=0
254            do m=1,nopem
255               if((nodedesiinv(nopesurf(m)).lt.0).or.
256     &            (nodedesiinv(nopesurf(m)).eq.1)) then
257                  nnodes=nnodes+1
258               endif
259            enddo
260!
261            if(nnodes.ge.nopedesi) then
262               do m=1,nopem
263                  if(nodedesiinv(nopesurf(m)).lt.0) then
264                     nodedesiinv(nopesurf(m))=1
265                  endif
266               enddo
267      else
268               do m=1,nopem
269                  if(nodedesiinv(nopesurf(m)).eq.-2) then
270                     nodedesiinv(nopesurf(m))=-1
271                  endif
272               enddo
273            endif
274
275            indexe=nodface(5,indexe)
276            if(indexe.eq.0) exit
277!
278         enddo
279      enddo
280!
281      kflag=1
282      call isortii(nodedesi,iaux,ndesi,kflag)
283!
284!------------------------------------------------------------------------
285!     if node i in nodedesiinv(i) is -2 or -1 --> delete node i from
286!     set of designvariables
287!------------------------------------------------------------------------
288!
289      do i=1,nk
290         if(nodedesiinv(i).eq.-2) then
291!
292            write(*,*) '*WARNING in getdesiinfo:'
293            write(*,*) '          node ',i,' is removed'
294            write(*,*) '          from the set of design'
295            write(*,*) '          variables as it is an'
296            write(*,*) '          internal node'
297            write(40,*) i
298            iwrite=1
299!
300         elseif(nodedesiinv(i).eq.-1) then
301!
302            write(*,*) '*WARNING in getdesiinfo:'
303            write(*,*) '          node ',i,' is removed'
304            write(*,*) '          from the set of design'
305            write(*,*) '          variables as not sufficient '
306            write(*,*) '          other variables are on the  '
307            write(*,*) '          surrounding element faces  '
308            write(40,*) i
309            iwrite=1
310!
311            nodedesiinv(i)=0
312            call nident(nodedesi,i,ndesi,id)
313            do k=id+1,ndesi
314               nodedesi(k-1)=nodedesi(k)
315            enddo
316            ndesi=ndesi-1
317         endif
318      enddo
319!
320!------------------------------------------------------------------------
321!     write remaining design variables in the dat file
322!------------------------------------------------------------------------
323!
324      write(5,*)'NUMBER OF DESIGN VARIABLES TAKEN INTO ACCOUNT '
325      write(5,*)'FOR RANDOM FIELD ASSESSMENT:'
326      write(5,*) ndesi
327      write(5,*)''
328      write(5,*)'NODE NUMBERS OF DESIGN VARIABLES: '
329      do i=1,ndesi
330         write(5,*) nodedesi(i),','
331      enddo
332!
333      if(iwrite.eq.1) then
334        write(*,*) '*INFO in getdesiinfo3d_robust:'
335        write(*,*) '      rejected design nodes are stored in'
336        write(*,*) '      file WarnNodeDesignReject.nam'
337        write(*,*) '      This file can be loaded into'
338        write(*,*) '      an active cgx-session by typing'
339        write(*,*)
340     &       '      read WarnNodeDesignReject.nam inp'
341        write(*,*)
342        close(40)
343      else
344        close(40,status='delete')
345      endif
346!
347      return
348      end
349
350