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

Go to the source code of this file.

Functions/Subroutines

subroutine map3dto1d2d (yn, ipkon, inum, kon, lakon, nfield, nk, ne, cflag, co, vold, force, mi)
 

Function/Subroutine Documentation

◆ map3dto1d2d()

subroutine map3dto1d2d ( real*8, dimension(nfield,*)  yn,
integer, dimension(*)  ipkon,
integer, dimension(*)  inum,
integer, dimension(*)  kon,
character*8, dimension(*)  lakon,
integer  nfield,
integer  nk,
integer  ne,
character*1  cflag,
real*8, dimension(3,*)  co,
real*8, dimension(0:mi(2),*)  vold,
logical  force,
integer, dimension(*)  mi 
)
21 !
22 ! interpolates 3d field nodal values to 1d/2d nodal locations
23 !
24 ! the number of internal state variables is limited to 999
25 ! (cfr. array field)
26 !
27  implicit none
28 !
29  logical force,quadratic
30 !
31  character*1 cflag
32  character*8 lakon(*),lakonl
33 !
34  integer ipkon(*),inum(*),kon(*),ne,indexe,nfield,nk,i,j,k,l,m,
35  & node3(8,3),node6(3,6),node8(3,8),node2d,node3d,indexe2d,ne1d2d,
36  & node3m(8,3),node(8),m1,m2,nodea,nodeb,nodec,iflag,mi(*),jmax,
37  & jinc,nope
38 !
39  real*8 yn(nfield,*),cg(3),p(3),pcg(3),t(3),xl(3,8),shp(7,8),
40  & xsj(3),e1(3),e2(3),e3(3),s(6),dd,xi,et,co(3,*),xs(3,7),
41  & vold(0:mi(2),*),ratioe(3)
42 !
43  include "gauss.f"
44 !
45  data node3 /1,4,8,5,12,20,16,17,9,11,15,13,
46  & 0,0,0,0,2,3,7,6,10,19,14,18/
47  data node3m /1,5,8,4,17,16,20,12,
48  & 0,0,0,0,0,0,0,0,
49  & 3,7,6,2,19,14,18,10/
50  data node6 /1,13,4,2,14,5,3,15,6,7,0,10,8,0,11,9,0,12/
51  data node8 /1,17,5,2,18,6,3,19,7,4,20,8,9,0,13,10,0,14,
52  & 11,0,15,12,0,16/
53  data ratioe /0.16666666666667d0,0.66666666666666d0,
54  & 0.16666666666667d0/
55  data iflag /2/
56 !
57 ! removing any results in 1d/2d nodes
58 !
59  ne1d2d=0
60 !
61  do i=1,ne
62 !
63  if(ipkon(i).lt.0) cycle
64  lakonl=lakon(i)
65  if((lakonl(7:7).eq.' ').or.(lakonl(7:7).eq.'I').or.
66  & (lakonl(1:1).ne.'C')) cycle
67  ne1d2d=1
68  indexe=ipkon(i)
69 !
70  if((lakonl(4:5).eq.'15').or.(lakonl(4:4).eq.'6')) then
71  if(lakonl(4:5).eq.'15') then
72  indexe2d=indexe+15
73  jmax=6
74  else
75  indexe2d=indexe+6
76  jmax=3
77  endif
78  do j=1,jmax
79  node2d=kon(indexe2d+j)
80  inum(node2d)=0
81  do k=1,nfield
82  yn(k,node2d)=0.d0
83  enddo
84  enddo
85  elseif(lakonl(7:7).eq.'B') then
86  if(lakonl(4:5).eq.'8I') then
87  indexe2d=indexe+11
88  jmax=2
89  elseif(lakonl(4:5).eq.'8R') then
90  indexe2d=indexe+8
91  jmax=2
92  elseif(lakonl(4:5).eq.'20') then
93  indexe2d=indexe+20
94  jmax=3
95  endif
96  do j=1,jmax
97  node2d=kon(indexe2d+j)
98  inum(node2d)=0
99  do k=1,nfield
100  yn(k,node2d)=0.d0
101  enddo
102  enddo
103  else
104  if(lakonl(4:5).eq.'8I') then
105  indexe2d=indexe+11
106  jmax=4
107  elseif((lakonl(4:5).eq.'8R').or.(lakonl(4:5).eq.'8 ')) then
108  indexe2d=indexe+8
109  jmax=4
110  elseif(lakonl(4:5).eq.'20') then
111  indexe2d=indexe+20
112  jmax=8
113  endif
114  do j=1,jmax
115  node2d=kon(indexe2d+j)
116  inum(node2d)=0
117  do k=1,nfield
118  yn(k,node2d)=0.d0
119  enddo
120  enddo
121  endif
122 !
123 ! inactivating the 3d expansion nodes of 1d/2d elements
124 ! in case forces are mapped this field is used to ensure
125 ! that the forces in the 3d-nodes are mapped only once onto the
126 ! 2d-nodes
127 !
128  do j=1,indexe2d-indexe
129  inum(kon(indexe+j))=0
130  enddo
131 !
132  enddo
133 !
134 ! if no 1d/2d elements return
135 !
136  if(ne1d2d.eq.0) return
137 !
138 ! interpolation of 3d results on 1d/2d nodes
139 !
140  do i=1,ne
141 !
142  if(ipkon(i).lt.0) cycle
143  lakonl=lakon(i)
144  if((lakonl(7:7).eq.' ').or.(lakonl(7:7).eq.'I').or.
145  & (lakonl(1:1).ne.'C')) cycle
146  indexe=ipkon(i)
147 !
148 ! check whether linear or quadratic element
149 !
150  if((lakonl(4:4).eq.'6').or.(lakonl(4:4).eq.'8')) then
151  quadratic=.false.
152  else
153  quadratic=.true.
154  endif
155 !
156  if((lakonl(4:5).eq.'15').or.(lakonl(4:4).eq.'6')) then
157  if(lakonl(4:5).eq.'15') then
158  indexe2d=indexe+15
159  jmax=6
160  else
161  indexe2d=indexe+6
162  jmax=3
163  endif
164  do j=1,jmax
165  node2d=kon(indexe2d+j)
166  inum(node2d)=inum(node2d)-1
167  if(.not.force) then
168 !
169 ! taking the mean across the thickness
170 !
171  if((j.le.3).and.(quadratic)) then
172 !
173 ! end nodes: weights 1/6,2/3 and 1/6
174 !
175  do l=1,3
176  node3d=kon(indexe+node6(l,j))
177  do k=1,nfield
178  yn(k,node2d)=yn(k,node2d)+
179  & yn(k,node3d)*ratioe(l)
180  enddo
181  enddo
182  else
183 !
184 ! middle nodes: weights 1/2,1/2
185 !
186  do l=1,3,2
187  node3d=kon(indexe+node6(l,j))
188  do k=1,nfield
189  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)/2.d0
190  enddo
191  enddo
192  endif
193  else
194 !
195 ! forces must be summed
196 !
197 ! the contribution of each 3d-node should only be taken
198 ! once. inum(node3d) is used as marker.
199 !
200  inum(node2d)=-1
201 !
202  if((j.le.3).and.(quadratic)) then
203 !
204 ! end nodes of quadratic 2d-elements
205 !
206  do l=1,3
207  node3d=kon(indexe+node6(l,j))
208  if(inum(node3d).ne.0) cycle
209  inum(node3d)=1
210  do k=1,nfield
211  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)
212  enddo
213  enddo
214  else
215 !
216 ! middle nodes of quadratic 2d-elements
217 ! or end nodes of linear 2d-elements
218 !
219  do l=1,3,2
220  node3d=kon(indexe+node6(l,j))
221  if(inum(node3d).ne.0) cycle
222  inum(node3d)=1
223  do k=1,nfield
224  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)
225  enddo
226  enddo
227  endif
228  endif
229  enddo
230  elseif(lakonl(7:7).eq.'B') then
231  if(lakonl(4:5).eq.'8I') then
232  indexe2d=indexe+11
233  jmax=2
234  jinc=1
235  nope=4
236  elseif(lakonl(4:5).eq.'8R') then
237  indexe2d=indexe+8
238  jmax=2
239  jinc=1
240  nope=4
241  elseif(lakonl(4:5).eq.'20') then
242  indexe2d=indexe+20
243  jmax=3
244  jinc=2
245  nope=8
246  endif
247  if(cflag.ne.'M') then
248 !
249 ! mean values for beam elements
250 !
251  do j=1,jmax
252  node2d=kon(indexe2d+j)
253  if(.not.force) then
254 !
255 ! mean value of vertex values
256 !
257  do l=1,4
258  inum(node2d)=inum(node2d)-1
259  if(quadratic) then
260  node3d=kon(indexe+node3(l,j))
261  else
262  node3d=kon(indexe+node3(l,2*j-1))
263  endif
264  do k=1,nfield
265  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)
266  enddo
267  enddo
268  else
269 !
270 ! forces must be summed across the section
271 !
272 ! the contribution of each 3d-node should only be taken
273 ! once. inum(node3d) is used as marker.
274 !
275  inum(node2d)=-1
276 !
277  if((j.ne.2).and.(quadratic)) then
278 !
279 ! end nodes of quadratic beam elements
280 !
281  do l=1,8
282  node3d=kon(indexe+node3(l,j))
283  if(inum(node3d).ne.0) cycle
284  inum(node3d)=1
285  do k=1,nfield
286  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)
287  enddo
288  enddo
289  else
290 !
291 ! middle nodes of quadratic beam elements or
292 ! end nodes of linear beam elements
293 !
294  do l=1,4
295  if(quadratic) then
296  node3d=kon(indexe+node3(l,j))
297  else
298  node3d=kon(indexe+node3(l,2*j-1))
299  endif
300  if(inum(node3d).ne.0) cycle
301  inum(node3d)=1
302  do k=1,nfield
303  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)
304  enddo
305  enddo
306  endif
307  endif
308  enddo
309  else
310 !
311 ! section forces for beam elements
312 !
313  do j=1,jmax,jinc
314  node2d=kon(indexe2d+j)
315  inum(node2d)=inum(node2d)-1
316 !
317 ! coordinates of the nodes belonging to the section
318 !
319  do l=1,nope
320  if(quadratic) then
321  node(l)=kon(indexe+node3m(l,j))
322  else
323  node(l)=kon(indexe+node3m(l,2*j-1))
324  endif
325  do m=1,3
326  xl(m,l)=co(m,node(l))+vold(m,node(l))
327  enddo
328  enddo
329 !
330 ! center of gravity and unit vectors 1 and 2
331 !
332  do m=1,3
333  cg(m)=(xl(m,1)+xl(m,2)+xl(m,3)+xl(m,4))/4.d0
334  if(j.eq.1) then
335  e1(m)=(xl(m,1)+xl(m,4)-xl(m,2)-xl(m,3))
336  else
337  e1(m)=(xl(m,2)+xl(m,3)-xl(m,1)-xl(m,4))
338  endif
339  e2(m)=(xl(m,3)+xl(m,4)-xl(m,1)-xl(m,2))
340  enddo
341 !
342 ! normalizing e1
343 !
344  dd=dsqrt(e1(1)*e1(1)+e1(2)*e1(2)+e1(3)*e1(3))
345  do m=1,3
346  e1(m)=e1(m)/dd
347  enddo
348 !
349 ! making sure that e2 is orthogonal to e1
350 !
351  dd=e1(1)*e2(1)+e1(2)*e2(2)+e1(3)*e2(3)
352  do m=1,3
353  e2(m)=e2(m)-dd*e1(m)
354  enddo
355 !
356 ! normalizing e2
357 !
358  dd=dsqrt(e2(1)*e2(1)+e2(2)*e2(2)+e2(3)*e2(3))
359  do m=1,3
360  e2(m)=e2(m)/dd
361  enddo
362 !
363 ! e3 = e1 x e2 for j=3, e3 = e2 x e1 for j=1
364 !
365  if(j.eq.1) then
366  e3(1)=e2(2)*e1(3)-e1(2)*e2(3)
367  e3(2)=e2(3)*e1(1)-e1(3)*e2(1)
368  e3(3)=e2(1)*e1(2)-e1(1)*e2(2)
369  else
370  e3(1)=e1(2)*e2(3)-e2(2)*e1(3)
371  e3(2)=e1(3)*e2(1)-e2(3)*e1(1)
372  e3(3)=e1(1)*e2(2)-e2(1)*e1(2)
373  endif
374 !
375 ! loop over the integration points (2x2)
376 !
377  do l=1,4
378  xi=gauss2d2(1,l)
379  et=gauss2d2(2,l)
380  if(quadratic) then
381  call shape8q(xi,et,xl,xsj,xs,shp,iflag)
382  else
383  call shape4q(xi,et,xl,xsj,xs,shp,iflag)
384  endif
385 !
386 ! local stress tensor
387 !
388  do m1=1,6
389  s(m1)=0.d0
390  do m2=1,nope
391  s(m1)=s(m1)+shp(4,m2)*yn(m1,node(m2))
392  enddo
393  enddo
394 !
395 ! local coordinates
396 !
397  do m1=1,3
398  p(m1)=0.d0
399  do m2=1,nope
400  p(m1)=p(m1)+shp(4,m2)*xl(m1,m2)
401  enddo
402  pcg(m1)=p(m1)-cg(m1)
403  enddo
404 !
405 ! local stress vector on section
406 !
407  t(1)=s(1)*xsj(1)+s(4)*xsj(2)+s(5)*xsj(3)
408  t(2)=s(4)*xsj(1)+s(2)*xsj(2)+s(6)*xsj(3)
409  t(3)=s(5)*xsj(1)+s(6)*xsj(2)+s(3)*xsj(3)
410 !
411 ! section forces
412 !
413  yn(1,node2d)=yn(1,node2d)+
414  & (e1(1)*t(1)+e1(2)*t(2)+e1(3)*t(3))
415  yn(2,node2d)=yn(2,node2d)+
416  & (e2(1)*t(1)+e2(2)*t(2)+e2(3)*t(3))
417  yn(3,node2d)=yn(3,node2d)+
418  & (e3(1)*t(1)+e3(2)*t(2)+e3(3)*t(3))
419 !
420 ! section moments
421 !
422 ! about beam axis
423 !
424  yn(4,node2d)=yn(4,node2d)+
425  & (e3(1)*pcg(2)*t(3)+e3(2)*pcg(3)*t(1)+
426  & e3(3)*pcg(1)*t(2)-e3(3)*pcg(2)*t(1)-
427  & e3(1)*pcg(3)*t(2)-e3(2)*pcg(1)*t(3))
428 !
429 ! about 2-direction
430 !
431  yn(5,node2d)=yn(5,node2d)+
432  & (e2(1)*pcg(2)*t(3)+e2(2)*pcg(3)*t(1)+
433  & e2(3)*pcg(1)*t(2)-e2(3)*pcg(2)*t(1)-
434  & e2(1)*pcg(3)*t(2)-e2(2)*pcg(1)*t(3))
435 !
436 ! about 1-direction
437 !
438  yn(6,node2d)=yn(6,node2d)+
439  & (e1(1)*pcg(2)*t(3)+e1(2)*pcg(3)*t(1)+
440  & e1(3)*pcg(1)*t(2)-e1(3)*pcg(2)*t(1)-
441  & e1(1)*pcg(3)*t(2)-e1(2)*pcg(1)*t(3))
442 !
443 ! components 5 and 6 are switched in the frd
444 ! format, so the final order is beam axis,
445 ! 1-direction and 2-direction, or s12, s23 and s31
446 !
447  enddo
448  enddo
449 !
450  endif
451  else
452  if(lakonl(4:5).eq.'8I') then
453  indexe2d=indexe+11
454  jmax=4
455  elseif((lakonl(4:5).eq.'8R').or.(lakonl(4:5).eq.'8 ')) then
456  indexe2d=indexe+8
457  jmax=4
458  elseif(lakonl(4:5).eq.'20') then
459  indexe2d=indexe+20
460  jmax=8
461  endif
462  do j=1,jmax
463  node2d=kon(indexe2d+j)
464  inum(node2d)=inum(node2d)-1
465  if(.not.force) then
466 !
467 ! taking the mean across the thickness
468 !
469  if((j.le.4).and.(quadratic)) then
470 !
471 ! end nodes: weights 1/6,2/3 and 1/6
472 !
473  do l=1,3
474  node3d=kon(indexe+node8(l,j))
475  do k=1,nfield
476  yn(k,node2d)=yn(k,node2d)+
477  & yn(k,node3d)*ratioe(l)
478  enddo
479  enddo
480  else
481 !
482 ! middle nodes: weights 1/2,1/2
483 !
484  do l=1,3,2
485  node3d=kon(indexe+node8(l,j))
486  do k=1,nfield
487  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)/2.d0
488  enddo
489  enddo
490  endif
491  else
492 !
493 ! forces must be summed
494 !
495 ! the contribution of each 3d-node should only be taken
496 ! once. inum(node3d) is used as marker.
497 !
498  inum(node2d)=-1
499 !
500  if((j.le.4).and.(quadratic)) then
501 !
502 ! end nodes of quadratic 2d-elements
503 !
504  do l=1,3
505  node3d=kon(indexe+node8(l,j))
506  if(inum(node3d).ne.0) cycle
507  inum(node3d)=1
508  do k=1,nfield
509  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)
510  enddo
511  enddo
512  else
513 !
514 ! middle nodes of quadratic 2d-elements
515 ! or end nodes of linear 2d-elements
516 !
517  do l=1,3,2
518  node3d=kon(indexe+node8(l,j))
519  if(inum(node3d).ne.0) cycle
520  inum(node3d)=1
521  do k=1,nfield
522  yn(k,node2d)=yn(k,node2d)+yn(k,node3d)
523  enddo
524  enddo
525  endif
526  endif
527  enddo
528  endif
529 !
530  enddo
531 !
532 ! taking the mean of nodal contributions coming from different
533 ! elements having the node in common
534 ! restoring inum for the 3d-nodes to zero
535 !
536  do i=1,ne
537 !
538  if(ipkon(i).lt.0) cycle
539  lakonl=lakon(i)
540  if((lakonl(7:7).eq.' ').or.(lakonl(7:7).eq.'I').or.
541  & (lakonl(1:1).ne.'C')) cycle
542  indexe=ipkon(i)
543 !
544  if((lakonl(4:5).eq.'15').or.(lakonl(4:4).eq.'6')) then
545  if(lakonl(4:5).eq.'15') then
546  indexe2d=indexe+15
547  jmax=6
548  else
549  indexe2d=indexe+6
550  jmax=3
551  endif
552  do j=1,jmax
553  node2d=kon(indexe2d+j)
554  if(inum(node2d).lt.0) then
555  inum(node2d)=-inum(node2d)
556  do k=1,nfield
557  yn(k,node2d)=yn(k,node2d)/inum(node2d)
558  enddo
559  endif
560  enddo
561  elseif(lakonl(7:7).eq.'B') then
562  if(lakonl(4:5).eq.'8I') then
563  indexe2d=indexe+11
564  jmax=2
565  elseif(lakonl(4:5).eq.'8R') then
566  indexe2d=indexe+8
567  jmax=2
568  elseif(lakonl(4:5).eq.'20') then
569  indexe2d=indexe+20
570  jmax=3
571  endif
572  do j=1,jmax
573  node2d=kon(indexe2d+j)
574  if(inum(node2d).lt.0) then
575  inum(node2d)=-inum(node2d)
576  do k=1,nfield
577  yn(k,node2d)=yn(k,node2d)/inum(node2d)
578  enddo
579  endif
580  enddo
581  else
582  if(lakonl(4:5).eq.'8I') then
583  indexe2d=indexe+11
584  jmax=4
585  elseif((lakonl(4:5).eq.'8R').or.(lakonl(4:5).eq.'8 ')) then
586  indexe2d=indexe+8
587  jmax=4
588  elseif(lakonl(4:5).eq.'20') then
589  indexe2d=indexe+20
590  jmax=8
591  endif
592  do j=1,jmax
593  node2d=kon(indexe2d+j)
594  if(inum(node2d).lt.0) then
595  inum(node2d)=-inum(node2d)
596  do k=1,nfield
597  yn(k,node2d)=yn(k,node2d)/inum(node2d)
598  enddo
599  endif
600  enddo
601  endif
602 !
603 ! inactivating the 3d expansion nodes of 1d/2d elements
604 ! in case forces are mapped this field is used to ensure
605 ! that the forces in the 3d-nodes are mapped only once onto the
606 ! 2d-nodes
607 !
608  do j=1,indexe2d-indexe
609  inum(kon(indexe+j))=0
610  enddo
611 !
612  enddo
613 !
614 ! beam section forces in the middle nodes
615 !
616  do i=1,ne
617 !
618  if(ipkon(i).lt.0) cycle
619  lakonl=lakon(i)
620  if((lakonl(7:7).eq.' ').or.(lakonl(7:7).eq.'I').or.
621  & (lakonl(1:1).ne.'C')) cycle
622  indexe=ipkon(i)
623 !
624 ! not relevant for linear elements
625 !
626  if((lakonl(4:4).eq.'6').or.(lakonl(4:4).eq.'8')) cycle
627 !
628  if(lakonl(7:7).eq.'B') then
629  indexe2d=indexe+20
630  if(cflag.eq.'M') then
631 !
632 ! section forces in the middle node are the mean
633 ! of those in the end nodes
634 !
635  nodea=kon(indexe2d+1)
636  nodeb=kon(indexe2d+2)
637  nodec=kon(indexe2d+3)
638  inum(nodeb)=1
639  do j=1,6
640  yn(j,nodeb)=yn(j,nodeb)+(yn(j,nodea)+yn(j,nodec))/2.d0
641  enddo
642 !
643  endif
644  endif
645 !
646  enddo
647 !
648  return
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)