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

Go to the source code of this file.

Functions/Subroutines

subroutine mafilldm (co, nk, kon, ipkon, lakon, ne, nodeboun, ndirboun, xboun, nboun, ipompc, nodempc, coefmpc, nmpc, nodeforc, ndirforc, xforc, nforc, nelemload, sideload, xload, nload, xbody, ipobody, nbody, cgr, ad, au, nactdof, icol, jq, irow, neq, nzl, nmethod, ikmpc, ilmpc, ikboun, ilboun, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, prestr, iprestr, vold, iperturb, sti, nzs, stx, adb, aub, iexpl, plicon, nplicon, plkcon, nplkcon, xstiff, npmat_, dtime, matname, mi, ncmat_, ttime, time, istep, iinc, ibody, clearini, mortar, springarea, pslavsurf, pmastsurf, reltime, nasym)
 

Function/Subroutine Documentation

◆ mafilldm()

subroutine mafilldm ( real*8, dimension(3,*)  co,
integer  nk,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
character*8, dimension(*)  lakon,
integer  ne,
integer, dimension(*)  nodeboun,
integer, dimension(*)  ndirboun,
real*8, dimension(*)  xboun,
integer  nboun,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer, dimension(2,*)  nodeforc,
integer, dimension(*)  ndirforc,
real*8, dimension(*)  xforc,
integer  nforc,
integer, dimension(2,*)  nelemload,
character*20, dimension(*)  sideload,
real*8, dimension(2,*)  xload,
integer  nload,
real*8, dimension(7,*)  xbody,
integer, dimension(2,*)  ipobody,
integer  nbody,
real*8, dimension(4,*)  cgr,
real*8, dimension(*)  ad,
real*8, dimension(*)  au,
integer, dimension(0:mi(2),*)  nactdof,
integer, dimension(*)  icol,
integer, dimension(*)  jq,
integer, dimension(*)  irow,
integer, dimension(2)  neq,
integer  nzl,
integer  nmethod,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
integer, dimension(*)  ikboun,
integer, dimension(*)  ilboun,
real*8, dimension(0:ncmat_,ntmat_,*)  elcon,
integer, dimension(2,*)  nelcon,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
real*8, dimension(0:6,ntmat_,*)  alcon,
integer, dimension(2,*)  nalcon,
real*8, dimension(*)  alzero,
integer, dimension(mi(3),*)  ielmat,
integer, dimension(mi(3),*)  ielorien,
integer  norien,
real*8, dimension(7,*)  orab,
integer  ntmat_,
real*8, dimension(*)  t0,
real*8, dimension(*)  t1,
integer  ithermal,
real*8, dimension(6,mi(1),*)  prestr,
integer  iprestr,
real*8, dimension(0:mi(2),*)  vold,
integer  iperturb,
real*8, dimension(6,mi(1),*)  sti,
integer, dimension(3)  nzs,
real*8, dimension(6,mi(1),*)  stx,
real*8, dimension(*)  adb,
real*8, dimension(*)  aub,
integer  iexpl,
real*8, dimension(0:2*npmat_,ntmat_,*)  plicon,
integer, dimension(0:ntmat_,*)  nplicon,
real*8, dimension(0:2*npmat_,ntmat_,*)  plkcon,
integer, dimension(0:ntmat_,*)  nplkcon,
real*8, dimension(27,mi(1),*)  xstiff,
integer  npmat_,
real*8  dtime,
character*80, dimension(*)  matname,
integer, dimension(*)  mi,
integer  ncmat_,
real*8  ttime,
real*8  time,
integer  istep,
integer  iinc,
integer, dimension(3,*)  ibody,
real*8, dimension(3,9,*)  clearini,
integer  mortar,
real*8, dimension(2,*)  springarea,
real*8, dimension(3,*)  pslavsurf,
real*8, dimension(6,*)  pmastsurf,
real*8  reltime,
integer  nasym 
)
31 !
32 ! filling the damping matrix in spare matrix format (sm)
33 !
34  implicit none
35 !
36  character*8 lakon(*)
37  character*20 sideload(*)
38  character*80 matname(*)
39 !
40  integer kon(*),nodeboun(*),ndirboun(*),ipompc(*),nodempc(3,*),
41  & nodeforc(2,*),ndirforc(*),nelemload(2,*),icol(*),jq(*),ikmpc(*),
42  & ilmpc(*),ikboun(*),ilboun(*),mi(*),nactdof(0:mi(2),*),konl(20),
43  & irow(*),nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),
44  & ielorien(mi(3),*),ipkon(*),ipobody(2,*),nbody,ibody(3,*),
45  & nk,ne,nboun,nmpc,nforc,nload,neq(2),nzl,nmethod,icolumn,
46  & ithermal,iprestr,iperturb,nzs(3),i,j,k,l,m,idist,jj,iloc,
47  & ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2,idof1,idof2,
48  & mpc1,mpc2,index1,index2,node1,node2,ne0,igauss,imat,
49  & ntmat_,indexe,nope,norien,iexpl,i0,ncmat_,istep,iinc,
50  & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,jfaces,
51  & kode,mortar,nasym
52 !
53  real*8 co(3,*),xboun(*),coefmpc(*),xforc(*),xload(2,*),p1(3),
54  & p2(3),ad(*),au(*),bodyf(3),voldl(0:mi(2),26),clearini(3,9,*),
55  & t0(*),t1(*),prestr(6,mi(1),*),vold(0:mi(2),*),s(60,60),
56  & ff(60),
57  & sti(6,mi(1),*),sm(60,60),stx(6,mi(1),*),adb(*),aub(*),
58  & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),elas(21),
59  & alcon(0:6,ntmat_,*),alzero(*),orab(7,*),xbody(7,*),cgr(4,*),
60  & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
61  & xstiff(27,mi(1),*),om,value,dtime,ttime,time,elconloc(21),
62  & xl(3,26),springarea(2,*),pslavsurf(3,*),reltime,pmastsurf(6,*)
63 !
64  i0=0
65 !
66 ! determining nzl
67 !
68  nzl=0
69  do i=neq(2),1,-1
70  if(icol(i).gt.0) then
71  nzl=i
72  exit
73  endif
74  enddo
75 c!
76 c! initializing the matrices
77 c!
78 c do i=1,neq(2)
79 c ad(i)=0.d0
80 c enddo
81 c do i=1,nzs(2)
82 c au(i)=0.d0
83 c enddo
84 !
85  if((ithermal.le.1).or.(ithermal.eq.3)) then
86 !
87 ! mechanical analysis: loop over all elements
88 !
89  ne0=0
90  do i=1,ne
91 !
92  if(ipkon(i).lt.0) cycle
93  if(lakon(i)(1:2).eq.'ED') then
94  indexe=ipkon(i)
95  read(lakon(i)(8:8),'(i1)') nope
96  nope=nope+1
97 !
98  do j=1,nope
99  konl(j)=kon(indexe+j)
100  enddo
101 !
102  call e_damp(co,nk,konl,lakon(i),p1,p2,om,bodyf,nbody,s,
103  & sm,ff,i,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
104  & alzero,ielmat,ielorien,norien,orab,ntmat_,
105  & t0,t1,ithermal,vold,iperturb,
106  & nelemload,sideload,xload,nload,idist,sti,stx,
107  & iexpl,plicon,nplicon,plkcon,nplkcon,xstiff,npmat_,
108  & dtime,matname,mi(1),ncmat_,ttime,time,istep,iinc,
109  & nmethod)
110  elseif((lakon(i)(1:2).eq.'ES').and.
111  & (lakon(i)(7:7).eq.'C')) then
112  indexe=ipkon(i)
113  if(mortar.eq.0) then
114  read(lakon(i)(8:8),'(i1)') nope
115  nope=nope+1
116  konl(nope+1)=kon(indexe+nope+1)
117  elseif(mortar.eq.1) then
118  nope=kon(indexe)
119  endif
120  imat=ielmat(1,i)
121 !
122 ! computation of the coordinates of the local nodes
123 !
124  do k=1,nope
125  konl(k)=kon(indexe+k)
126  do j=1,3
127  xl(j,k)=co(j,konl(k))
128  voldl(j,k)=vold(j,konl(k))
129  enddo
130  enddo
131 !
132 ! initialisation of s
133 !
134  do k=1,3*nope
135  do j=1,3*nope
136  s(k,j)=0.d0
137  enddo
138  enddo
139 !
140  kode=nelcon(1,imat)
141 !
142 ! as soon as the first contact element is discovered ne0 is
143 ! determined and saved
144 !
145  if(ne0.eq.0) ne0=i-1
146  if(mortar.eq.0) then
147  call springdamp_n2f(xl,elas,voldl,s,imat,elcon,
148  & ncmat_,ntmat_,nope,iperturb,
149  & springarea(1,konl(nope+1)),nmethod,
150  & mi,reltime,nasym)
151  elseif(mortar.eq.1) then
152  iloc=kon(indexe+nope+1)
153  jfaces=kon(indexe+nope+2)
154  igauss=kon(indexe+nope+1)
155  call springdamp_f2f(xl,elas,voldl,s,imat,elcon,
156  & ncmat_,ntmat_,nope,lakon(i),iperturb,
157  & springarea(1,iloc),
158  & nmethod,mi,reltime,nasym,jfaces,igauss,pslavsurf,
159  & pmastsurf,clearini)
160  endif
161  else
162  cycle
163  endif
164 !
165  if(nasym.eq.0) then
166 !
167  do jj=1,3*nope
168 !
169  j=(jj-1)/3+1
170  k=jj-3*(j-1)
171 !
172  node1=kon(indexe+j)
173  jdof1=nactdof(k,node1)
174 !
175  do ll=jj,3*nope
176 !
177  l=(ll-1)/3+1
178  m=ll-3*(l-1)
179 !
180  node2=kon(indexe+l)
181  jdof2=nactdof(m,node2)
182 !
183 ! check whether one of the DOF belongs to a SPC or MPC
184 !
185  if((jdof1.gt.0).and.(jdof2.gt.0)) then
186  call add_sm_st(au,ad,jq,irow,jdof1,jdof2,
187  & s(jj,ll),jj,ll)
188  elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
189 !
190 ! idof1: genuine DOF
191 ! idof2: nominal DOF of the SPC/MPC
192 !
193  if(jdof1.le.0) then
194  idof1=jdof2
195 c idof2=(node1-1)*8+k
196  idof2=jdof1
197  else
198  idof1=jdof1
199 c idof2=(node2-1)*8+m
200  idof2=jdof2
201  endif
202  if(nmpc.gt.0) then
203 c call nident(ikmpc,idof2,nmpc,id)
204 c if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
205  if(idof2.ne.2*(idof2/2)) then
206 !
207 ! regular DOF / MPC
208 !
209 c id=ilmpc(id)
210  id=(-idof2+1)/2
211  ist=ipompc(id)
212  index=nodempc(3,ist)
213  if(index.eq.0) cycle
214  do
215  idof2=nactdof(nodempc(2,index),
216  & nodempc(1,index))
217  value=-coefmpc(index)*s(jj,ll)/
218  & coefmpc(ist)
219  if(idof1.eq.idof2) value=2.d0*value
220  if(idof2.gt.0) then
221  call add_sm_st(au,ad,jq,irow,idof1,
222  & idof2,value,i0,i0)
223  endif
224  index=nodempc(3,index)
225  if(index.eq.0) exit
226  enddo
227  cycle
228  endif
229  endif
230  else
231 c idof1=(node1-1)*8+k
232 c idof2=(node2-1)*8+m
233  idof1=jdof1
234  idof2=jdof2
235  mpc1=0
236  mpc2=0
237  if(nmpc.gt.0) then
238 c call nident(ikmpc,idof1,nmpc,id1)
239 c if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
240 c call nident(ikmpc,idof2,nmpc,id2)
241 c if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
242  if(idof1.ne.2*(idof1/2)) mpc1=1
243  if(idof2.ne.2*(idof2/2)) mpc2=1
244  endif
245  if((mpc1.eq.1).and.(mpc2.eq.1)) then
246 c id1=ilmpc(id1)
247 c id2=ilmpc(id2)
248  id1=(-idof1+1)/2
249  id2=(-idof2+1)/2
250  if(id1.eq.id2) then
251 !
252 ! MPC id1 / MPC id1
253 !
254  ist=ipompc(id1)
255  index1=nodempc(3,ist)
256  if(index1.eq.0) cycle
257  do
258  idof1=nactdof(nodempc(2,index1),
259  & nodempc(1,index1))
260  index2=index1
261  do
262  idof2=nactdof(nodempc(2,index2),
263  & nodempc(1,index2))
264  value=coefmpc(index1)*coefmpc(index2)*
265  & s(jj,ll)/coefmpc(ist)/coefmpc(ist)
266  if((idof1.gt.0).and.(idof2.gt.0)) then
267  call add_sm_st(au,ad,jq,irow,
268  & idof1,idof2,value,i0,i0)
269  endif
270 !
271  index2=nodempc(3,index2)
272  if(index2.eq.0) exit
273  enddo
274  index1=nodempc(3,index1)
275  if(index1.eq.0) exit
276  enddo
277  else
278 !
279 ! MPC id1 / MPC id2
280 !
281  ist1=ipompc(id1)
282  index1=nodempc(3,ist1)
283  if(index1.eq.0) cycle
284  do
285  idof1=nactdof(nodempc(2,index1),
286  & nodempc(1,index1))
287  ist2=ipompc(id2)
288  index2=nodempc(3,ist2)
289  if(index2.eq.0) then
290  index1=nodempc(3,index1)
291  if(index1.eq.0) then
292  exit
293  else
294  cycle
295  endif
296  endif
297  do
298  idof2=nactdof(nodempc(2,index2),
299  & nodempc(1,index2))
300  value=coefmpc(index1)*coefmpc(index2)*
301  & s(jj,ll)/coefmpc(ist1)/
302  & coefmpc(ist2)
303  if(idof1.eq.idof2) value=2.d0*value
304  if((idof1.gt.0).and.(idof2.gt.0)) then
305  call add_sm_st(au,ad,jq,irow,
306  & idof1,idof2,value,i0,i0)
307  endif
308 !
309  index2=nodempc(3,index2)
310  if(index2.eq.0) exit
311  enddo
312  index1=nodempc(3,index1)
313  if(index1.eq.0) exit
314  enddo
315  endif
316  endif
317  endif
318  enddo
319 !
320  enddo
321 !
322 ! asymmetrical calculations (due to friction or tangential damping...)
323 !
324  else
325 !
326  do jj=1,3*nope
327 !
328  j=(jj-1)/3+1
329  k=jj-3*(j-1)
330 !
331  node1=kon(indexe+j)
332  jdof1=nactdof(k,node1)
333 !
334  do ll=1,3*nope
335 !
336  l=(ll-1)/3+1
337  m=ll-3*(l-1)
338 !
339  node2=kon(indexe+l)
340  jdof2=nactdof(m,node2)
341 !
342 ! check whether one of the DOF belongs to a SPC or MPC
343 !
344  if((jdof1.gt.0).and.(jdof2.gt.0)) then
345  call add_sm_st_as(au,ad,jq,irow,jdof1,jdof2,
346  & s(jj,ll),jj,ll,nzs)
347  elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
348 !
349 ! idof1: genuine DOF
350 ! idof2: nominal DOF of the SPC/MPC
351 !
352  if(jdof1.le.0) then
353 c idof1=(node1-1)*8+k
354  idof1=jdof1
355  idof2=jdof2
356  if(nmpc.gt.0) then
357 c call nident(ikmpc,idof1,nmpc,id)
358 c if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
359  if(idof1.ne.2*(idof1/2)) then
360 !
361 ! regular DOF / MPC
362 !
363 c id=ilmpc(id)
364  id=(-idof1+1)/2
365  ist=ipompc(id)
366  index=nodempc(3,ist)
367  if(index.eq.0) cycle
368  do
369  idof1=nactdof(nodempc(2,index),
370  & nodempc(1,index))
371  value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
372  if(idof1.gt.0) then
373  call add_sm_st_as(au,ad,jq,irow,idof1,
374  & idof2,value,i0,i0,nzs)
375  endif
376  index=nodempc(3,index)
377  if(index.eq.0) exit
378  enddo
379  cycle
380  endif
381  endif
382  else
383  idof1=jdof1
384 c idof2=(node2-1)*8+m
385  idof2=jdof2
386  if(nmpc.gt.0) then
387 c call nident(ikmpc,idof2,nmpc,id)
388 c if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
389  if(idof2.ne.2*(idof2/2)) then
390 !
391 ! regular DOF / MPC
392 !
393 c id=ilmpc(id)
394  id=(-idof2+1)/2
395  ist=ipompc(id)
396  index=nodempc(3,ist)
397  if(index.eq.0) cycle
398  do
399  idof2=nactdof(nodempc(2,index),
400  & nodempc(1,index))
401  value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
402  if(idof2.gt.0) then
403  call add_sm_st_as(au,ad,jq,irow,idof1,
404  & idof2,value,i0,i0,nzs)
405  endif
406  index=nodempc(3,index)
407  if(index.eq.0) exit
408  enddo
409  cycle
410  endif
411  endif
412  endif
413  else
414 c idof1=(node1-1)*8+k
415 c idof2=(node2-1)*8+m
416  idof1=jdof1
417  idof2=jdof2
418  mpc1=0
419  mpc2=0
420  if(nmpc.gt.0) then
421 c call nident(ikmpc,idof1,nmpc,id1)
422 c if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
423 c call nident(ikmpc,idof2,nmpc,id2)
424 c if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
425  if(idof1.ne.2*(idof1/2)) mpc1=1
426  if(idof2.ne.2*(idof2/2)) mpc2=1
427  endif
428  if((mpc1.eq.1).and.(mpc2.eq.1)) then
429 c id1=ilmpc(id1)
430 c id2=ilmpc(id2)
431  id1=(-idof1+1)/2
432  id2=(-idof2+1)/2
433  if(id1.eq.id2) then
434 !
435 ! MPC id1 / MPC id1
436 !
437  ist1=ipompc(id1)
438  index1=nodempc(3,ist1)
439  if(index1.eq.0) cycle
440  do
441  idof1=nactdof(nodempc(2,index1),
442  & nodempc(1,index1))
443 c index2=index1
444  ist2=ipompc(id2)
445  index2=nodempc(3,ist2)
446  if(index2.eq.0) then
447  index1=nodempc(3,index1)
448  if(index1.eq.0) then
449  exit
450  else
451  cycle
452  endif
453  endif
454  do
455  idof2=nactdof(nodempc(2,index2),
456  & nodempc(1,index2))
457 c value=coefmpc(index1)*coefmpc(index2)*
458 c & s(jj,ll)/coefmpc(ist)/coefmpc(ist)
459  value=coefmpc(index1)*coefmpc(index2)*
460  & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
461  if((idof1.gt.0).and.(idof2.gt.0)) then
462  call add_sm_st_as(au,ad,jq,irow,
463  & idof1,idof2,value,i0,i0,nzs)
464  endif
465 !
466  index2=nodempc(3,index2)
467  if(index2.eq.0) exit
468  enddo
469  index1=nodempc(3,index1)
470  if(index1.eq.0) exit
471  enddo
472  else
473 !
474 ! MPC id1 / MPC id2
475 !
476  ist1=ipompc(id1)
477  index1=nodempc(3,ist1)
478  if(index1.eq.0) cycle
479  do
480  idof1=nactdof(nodempc(2,index1),
481  & nodempc(1,index1))
482  ist2=ipompc(id2)
483  index2=nodempc(3,ist2)
484  if(index2.eq.0) then
485  index1=nodempc(3,index1)
486  if(index1.eq.0) then
487  exit
488  else
489  cycle
490  endif
491  endif
492  do
493  idof2=nactdof(nodempc(2,index2),
494  & nodempc(1,index2))
495  value=coefmpc(index1)*coefmpc(index2)*
496  & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
497  if((idof1.gt.0).and.(idof2.gt.0)) then
498  call add_sm_st_as(au,ad,jq,irow,
499  & idof1,idof2,value,i0,i0,nzs)
500  endif
501 !
502  index2=nodempc(3,index2)
503  if(index2.eq.0) exit
504  enddo
505  index1=nodempc(3,index1)
506  if(index1.eq.0) exit
507  enddo
508  endif
509  endif
510  endif
511  enddo
512  enddo
513  endif
514  enddo
515 !
516  endif
517 !
518 c do i=1,neq(2)
519 c write(*,*) i,ad(i)
520 c enddo
521 c do i=1,nzs(2)
522 c write(*,*) i,au(i)
523 c enddo
524  return
springdamp_f2f
subroutine springdamp_f2f(xl, elas, voldl, s, imat, elcon, ncmat_, ntmat_, nope, lakonl, iperturb, springarea, nmethod, mi, reltime, nasym, jfaces, igauss, pslavsurf, pmastsurf, clearini)
Definition: springdamp_f2f.f:22
e_damp
subroutine e_damp(co, nk, konl, lakonl, p1, p2, omx, bodyfx, nbody, s, sm, ff, nelem, elcon, nelcon, rhcon, nrhcon, alcon, nalcon, alzero, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, vold, iperturb, nelemload, sideload, xload, nload, idist, sti, stx, iexpl, plicon, nplicon, plkcon, nplkcon, xstiff, npmat_, dtime, matname, mi, ncmat_, ttime, time, istep, iinc, nmethod)
Definition: e_damp.f:26
springdamp_n2f
subroutine springdamp_n2f(xl, elas, voldl, s, imat, elcon, ncmat_, ntmat_, nope, iperturb, springarea, nmethod, mi, reltime, nasym)
Definition: springdamp_n2f.f:22
idist
static ITG * idist
Definition: radflowload.c:39
add_sm_st
subroutine add_sm_st(au, ad, jq, irow, i, j, value, i0, i1)
Definition: add_sm_st.f:20
add_sm_st_as
subroutine add_sm_st_as(au, ad, jq, irow, i, j, value, i0, i1, nzs)
Definition: add_sm_st_as.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)