CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
formgradient.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine formgradient (istartdesi, ialdesi, ipkon, lakon, ipoface, ndesi, nodedesi, nodface, kon, co, dgdx, nobject, weightformgrad, nodedesiinv, noregion, objectset, dgdxglob, nk)
 

Function/Subroutine Documentation

◆ formgradient()

subroutine formgradient ( integer, dimension(*)  istartdesi,
integer, dimension(*)  ialdesi,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer, dimension(*)  ipoface,
integer  ndesi,
integer, dimension(*)  nodedesi,
integer, dimension(5,*)  nodface,
integer, dimension(*)  kon,
real*8, dimension(3,*)  co,
real*8, dimension(ndesi,nobject)  dgdx,
integer  nobject,
real*8, dimension(ndesi)  weightformgrad,
integer, dimension(*)  nodedesiinv,
integer  noregion,
character*81, dimension(4,*)  objectset,
real*8, dimension(2,nk,nobject)  dgdxglob,
integer  nk 
)
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(4,*)
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,iaux,ifour,kflag,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,noregion,iactnode
43 !
44  real*8 xi,et,weight,xl(3,9),xs(3,2),xsj(3),shp(7,9),
45  & co(3,*),xsjj,weightformgrad(ndesi),dgdx(ndesi,nobject),
46  & dgdxglob(2,nk,nobject)
47 !
48 ! flag for shape functions
49 !
50  data iflag /3/
51  data indexel /0/
52 !
53  save indexel
54 !
55  include "gauss.f"
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 ! ifacq indicates for each local node number to which
88 ! faces this node belongs and the internal number
89 ! within this face. e.g., local node number 1 belongs to
90 ! face 1 with internal number 4, face 3 with internal number
91 ! 1 and face 6 with internal number 2 (first line underneath)
92 !
93 ! ifacq, ifact and ifacw can be derived in a unique way from
94 ! ifaceq, ifacet and ifacew2
95 !
96  data ifacq /1,4,3,1,6,2,
97  &1,3,3,2,4,1,
98  &1,2,4,2,5,1,
99  &1,1,5,2,6,1,
100  &2,1,3,4,6,3,
101  &2,2,3,3,4,4,
102  &2,3,4,3,5,4,
103  &2,4,5,3,6,4,
104  &1,7,3,5,0,0,
105  &1,6,4,5,0,0,
106  &1,5,5,5,0,0,
107  &1,8,6,5,0,0,
108  &2,5,3,7,0,0,
109  &2,6,4,7,0,0,
110  &2,7,5,7,0,0,
111  &2,8,6,7,0,0,
112  &3,8,6,6,0,0,
113  &3,6,4,8,0,0,
114  &4,6,5,8,0,0,
115  &5,6,6,8,0,0/
116 !
117  data ifact /1,1,2,1,4,1,
118  &1,3,2,2,3,1,
119  &1,2,3,2,4,3,
120  &2,3,3,3,4,2,
121  &1,6,2,4,0,0,
122  &1,5,3,4,0,0,
123  &1,4,4,6,0,0,
124  &2,6,4,4,0,0,
125  &2,5,3,6,0,0,
126  &3,5,4,5,0,0/
127 !
128  data ifacw /1,1,3,1,5,2,
129  &1,3,3,2,4,1,
130  &1,2,4,2,5,1,
131  &2,1,3,4,5,3,
132  &2,2,3,3,4,4,
133  &2,3,4,3,5,4,
134  &1,6,3,5,0,0,
135  &1,5,4,5,0,0,
136  &1,4,5,5,0,0,
137  &2,4,3,7,0,0,
138  &2,5,4,7,0,0,
139  &2,6,5,7,0,0,
140  &3,8,5,6,0,0,
141  &3,6,4,8,0,0,
142  &4,6,5,8,0,0/
143 !
144 ! flag for shape functions
145 !
146  kflag=1
147  ifour=4
148  ithree=3
149 !
150 !
151 ! Loop over all designvariables
152 !
153  do idesvar=1,ndesi
154 !
155  node=nodedesi(idesvar)
156 !
157 ! Loop over all elements belonging to the designvariable
158 !
159  do j=istartdesi(idesvar),istartdesi(idesvar+1)-1
160  ielem=ialdesi(j)
161 
162  indexe=ipkon(ielem)
163 !
164 ! mint2d: # of integration points on the surface
165 ! nopes: # of nodes in the surface
166 ! nope: # of nodes in the element
167 !
168  if(lakon(ielem)(4:5).eq.'8R') then
169  mint2d=1
170  nopes=4
171  nope=8
172  nopedesi=3
173  elseif(lakon(ielem)(4:4).eq.'8') then
174  mint2d=4
175  nopes=4
176  nope=8
177  nopedesi=3
178  elseif(lakon(ielem)(4:5).eq.'20') then
179  mint2d=9
180  nopes=8
181  nope=20
182  nopedesi=5
183  elseif(lakon(ielem)(4:5).eq.'10') then
184  mint2d=3
185  nopes=6
186  nope=10
187  nopedesi=3
188  elseif(lakon(ielem)(4:4).eq.'4') then
189  mint2d=1
190  nopes=3
191  nope=4
192  nopedesi=3
193  elseif(lakon(ielem)(4:4).eq.'6') then
194  mint2d=1
195  nopes=3
196  nope=6
197  nopedesi=3
198  elseif(lakon(ielem)(4:5).eq.'15') then
199  mint2d=3
200  nopes=6
201  nope=15
202  nopedesi=3
203  else
204  exit
205  endif
206  if(noregion.eq.1) nopedesi=0
207 !
208 ! actual position of the nodes belonging to the
209 ! master surface
210 !
211  do l=1,nope
212  konl(l)=kon(indexe+l)
213  enddo
214 !
215  do i=1,nope
216  if(kon(indexe+i).eq.node) exit
217  enddo
218 !
219 ! Loop over all faces of the actual element
220 !
221  loop1: do m=1,3
222 !
223 ! hexahedral elements
224 !
225  if((nope.eq.20).or.(nope.eq.8)) then
226  k=ifacq(1,m,i)
227  m1=ifacq(2,m,i)
228  do l=1,nopes
229  nodes1(l)=konl(ifaceq(l,k))
230  nodes2(l)=nodes1(l)
231  enddo
232  call isortii(nodes1,iaux,ifour,kflag)
233 !
234 ! tetrahedral element
235 !
236  elseif((nope.eq.10).or.(nope.eq.4)) then
237  k=ifact(1,m,i)
238  m1=ifact(2,m,i)
239  do l=1,nopes
240  nodes1(l)=konl(ifacet(l,k))
241  nodes2(l)=nodes1(l)
242  enddo
243  call isortii(nodes1,iaux,ithree,kflag)
244  nodes1(4)=0
245 !
246 ! wedge element
247 !
248  elseif(nope.eq.15) then
249  k=ifacw(1,m,i)
250  m1=ifacw(2,m,i)
251  if(k.le.2) then
252  nopes=6
253  mint3d=3
254  do l=1,nopes
255  nodes1(l)=konl(ifacew2(l,k))
256  nodes2(l)=nodes1(l)
257  enddo
258  call isortii(nodes1,iaux,ithree,kflag)
259  nodes1(4)=0
260  else
261  nopes=8
262  mint3d=4
263  do l=1,nopes
264  nodes1(l)=konl(ifacew2(l,k))
265  nodes2(l)=nodes1(l)
266  enddo
267  call isortii(nodes1,iaux,ithree,kflag)
268  endif
269  else
270  k=ifacw(1,m,i)
271  m1=ifacw(2,m,i)
272  if(k.le.2) then
273  nopes=3
274  mint3d=1
275  do l=1,nopes
276  nodes1(l)=konl(ifacew1(l,k))
277  nodes2(l)=nodes1(l)
278  enddo
279  call isortii(nodes1,iaux,ithree,kflag)
280  nodes1(4)=0
281  else
282  nopes=4
283  mint3d=2
284  do l=1,nopes
285  nodes1(l)=konl(ifacew1(l,k))
286  nodes2(l)=nodes1(l)
287  enddo
288  call isortii(nodes1,iaux,ithree,kflag)
289  endif
290  endif
291  if(k.eq.0) exit
292 !
293 ! Find external element face
294 !
295  indexs=ipoface(nodes1(1))
296  do
297  if(indexs.eq.0) cycle loop1
298  if((nodface(1,indexs).eq.nodes1(2)).and.
299  & (nodface(2,indexs).eq.nodes1(3))) then
300 !
301 ! Check if sufficient designnodes on that surface
302 !
303  nnodes=0
304  do n=1,nopes
305  if(nodedesiinv(nodes1(n)).eq.1) then
306  nnodes=nnodes+1
307  endif
308  enddo
309 !
310 ! Assign coordinates to nodes of surface
311 !
312  if(nnodes.ge.nopedesi) then
313  if((nope.eq.20).or.(nope.eq.8)) then
314  do l=1,nopes
315  do n=1,3
316  ifacel=ifaceq(l,nodface(4,indexs))
317  xl(n,l)=co(n,konl(ifacel))
318  enddo
319  enddo
320  elseif((nope.eq.10).or.(nope.eq.4)) then
321  do l=1,nopes
322  do n=1,3
323  ifacel=ifacet(l,nodface(4,indexs))
324  xl(n,l)=co(n,konl(ifacel))
325  enddo
326  enddo
327  elseif(nope.eq.15) then
328  do l=1,nopes
329  do n=1,3
330  ifacel=ifacew2(l,nodface(4,indexs))
331  xl(n,l)=co(n,konl(ifacel))
332  enddo
333  enddo
334  elseif(nope.eq.6) then
335  do l=1,nopes
336  do n=1,3
337  ifacel=ifacew1(l,nodface(4,indexs))
338  xl(n,l)=co(n,konl(ifacel))
339  enddo
340  enddo
341  endif
342 !
343 ! Determine weighting factors
344 !
345  do kk=1,mint2d
346  if(lakon(ielem)(4:5).eq.'20') then
347  xi=gauss2d3(1,kk)
348  et=gauss2d3(2,kk)
349  weight=weight2d3(kk)
350  call shape8q(xi,et,xl,xsj,
351  & xs,shp,iflag)
352  elseif(lakon(ielem)(4:5).eq.'10') then
353  xi=gauss2d5(1,kk)
354  et=gauss2d5(2,kk)
355  weight=weight2d5(kk)
356  call shape6tri(xi,et,xl,xsj,
357  & xs,shp,iflag)
358  elseif(lakon(ielem)(4:4).eq.'4') then
359  xi=gauss2d4(1,kk)
360  et=gauss2d4(2,kk)
361  weight=weight2d4(kk)
362  call shape3tri(xi,et,xl,xsj,
363  & xs,shp,iflag)
364  elseif(lakon(ielem)(4:5).eq.'15') then
365  if(k.le.2) then
366  xi=gauss2d5(1,kk)
367  et=gauss2d5(2,kk)
368  weight=weight2d5(kk)
369  call shape6tri(xi,et,xl,xsj,
370  & xs,shp,iflag)
371  else
372  xi=gauss2d3(1,kk)
373  et=gauss2d3(2,kk)
374  weight=weight2d3(kk)
375  call shape8q(xi,et,xl,xsj,
376  & xs,shp,iflag)
377  endif
378  elseif(lakon(ielem)(4:4).eq.'8') then
379  xi=gauss2d2(1,kk)
380  et=gauss2d2(2,kk)
381  weight=weight2d2(kk)
382  call shape4q(xi,et,xl,xsj,
383  & xs,shp,iflag)
384  elseif(lakon(ielem)(4:4).eq.'6') then
385  if(k.le.2) then
386  xi=gauss2d4(1,kk)
387  et=gauss2d4(2,kk)
388  weight=weight2d4(kk)
389  call shape3tri(xi,et,xl,xsj,
390  & xs,shp,iflag)
391  else
392  xi=gauss2d2(1,kk)
393  et=gauss2d2(2,kk)
394  weight=weight2d2(kk)
395  call shape4q(xi,et,xl,xsj,
396  & xs,shp,iflag)
397  endif
398  endif
399 !
400 ! Calculate Jacobian determinant
401 !
402  xsjj=dsqrt(xsj(1)**2+xsj(2)**2+
403  & xsj(3)**2)
404 !
405 ! Evaluate Shape functions for weightmatrix
406 !
407  weightformgrad(idesvar)=weightformgrad
408  & (idesvar)+weight*shp(4,m1)*xsjj
409  enddo
410  endif
411  exit
412  endif
413  indexs=nodface(5,indexs)
414  enddo ! find external element faces
415 !
416  enddo loop1 ! Loop over all faces of the actual element
417 !
418  enddo ! Loop over all elements belonging to the designvariable
419 !
420  enddo ! Loop over all designvariables
421 !
422 ! Scaling of sensitivities
423 !
424  do idesvar=1,ndesi
425  do iobject=1,nobject
426  if(objectset(1,iobject)(1:9).eq.'THICKNESS') cycle
427  if(weightformgrad(idesvar).gt.0.d0) then
428  dgdx(idesvar,iobject)=dgdx(idesvar,iobject)
429  & /weightformgrad(idesvar)
430  iactnode=nodedesi(idesvar)
431  dgdxglob(1,iactnode,iobject)=dgdx(idesvar,iobject)
432  else
433  iactnode=nodedesi(idesvar)
434  dgdxglob(1,iactnode,iobject)=dgdx(idesvar,iobject)
435  endif
436  enddo
437  enddo
438 !
439  return
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
isortii
subroutine isortii(ix, iy, n, kflag)
Definition: isortii.f:6
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)