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

Go to the source code of this file.

Functions/Subroutines

subroutine mafillsm (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, fext, 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_, mass, stiffness, buckling, rhsi, intscheme, physcon, shcon, nshcon, cocon, ncocon, ttime, time, istep, iinc, coriolis, ibody, xloadold, reltime, veold, springarea, nstate_, xstateini, xstate, thicke, integerglob, doubleglob, tieset, istartset, iendset, ialset, ntie, nasym, pslavsurf, pmastsurf, mortar, clearini, ielprop, prop, ne0, fnext, nea, neb, kscale, iponoel, inoel, network)
 

Function/Subroutine Documentation

◆ mafillsm()

subroutine mafillsm ( real*8, dimension(3,*), intent(in)  co,
integer, intent(in)  nk,
integer, dimension(*), intent(in)  kon,
integer, dimension(*), intent(in)  ipkon,
character*8, dimension(*), intent(in)  lakon,
integer, intent(in)  ne,
integer, dimension(*), intent(in)  nodeboun,
integer, dimension(*), intent(in)  ndirboun,
real*8, dimension(*), intent(in)  xboun,
integer, intent(in)  nboun,
integer, dimension(*), intent(in)  ipompc,
integer, dimension(3,*), intent(in)  nodempc,
real*8, dimension(*), intent(in)  coefmpc,
integer, intent(in)  nmpc,
integer, dimension(2,*), intent(in)  nodeforc,
integer, dimension(*), intent(in)  ndirforc,
real*8, dimension(*), intent(in)  xforc,
integer, intent(in)  nforc,
integer, dimension(2,*), intent(in)  nelemload,
character*20, dimension(*), intent(in)  sideload,
real*8, dimension(2,*), intent(inout)  xload,
integer, intent(in)  nload,
real*8, dimension(7,*), intent(in)  xbody,
integer, dimension(2,*), intent(in)  ipobody,
integer, intent(in)  nbody,
real*8, dimension(4,*), intent(inout)  cgr,
real*8, dimension(*), intent(inout)  ad,
real*8, dimension(*), intent(inout)  au,
real*8, dimension(*), intent(inout)  fext,
integer, dimension(0:mi(2),*), intent(in)  nactdof,
integer, dimension(*), intent(in)  icol,
integer, dimension(*), intent(in)  jq,
integer, dimension(*), intent(in)  irow,
integer, dimension(2), intent(in)  neq,
integer, intent(in)  nzl,
integer, intent(inout)  nmethod,
integer, dimension(*), intent(in)  ikmpc,
integer, dimension(*), intent(in)  ilmpc,
integer, dimension(*), intent(in)  ikboun,
integer, dimension(*), intent(in)  ilboun,
real*8, dimension(0:ncmat_,ntmat_,*), intent(in)  elcon,
integer, dimension(2,*), intent(in)  nelcon,
real*8, dimension(0:1,ntmat_,*), intent(in)  rhcon,
integer, dimension(*), intent(in)  nrhcon,
real*8, dimension(0:6,ntmat_,*), intent(in)  alcon,
integer, dimension(2,*), intent(in)  nalcon,
real*8, dimension(*), intent(in)  alzero,
integer, dimension(mi(3),*), intent(in)  ielmat,
integer, dimension(mi(3),*), intent(in)  ielorien,
integer, intent(in)  norien,
real*8, dimension(7,*), intent(in)  orab,
integer, intent(in)  ntmat_,
real*8, dimension(*), intent(in)  t0,
real*8, dimension(*), intent(in)  t1,
integer, dimension(2), intent(in)  ithermal,
real*8, dimension(6,mi(1),*), intent(in)  prestr,
integer, intent(in)  iprestr,
real*8, dimension(0:mi(2),*), intent(in)  vold,
integer, dimension(*), intent(in)  iperturb,
real*8, dimension(6,mi(1),*), intent(in)  sti,
integer, dimension(3), intent(in)  nzs,
real*8, dimension(6,mi(1),*), intent(in)  stx,
real*8, dimension(*), intent(inout)  adb,
real*8, dimension(*), intent(inout)  aub,
integer, intent(in)  iexpl,
real*8, dimension(0:2*npmat_,ntmat_,*), intent(in)  plicon,
integer, dimension(0:ntmat_,*), intent(in)  nplicon,
real*8, dimension(0:2*npmat_,ntmat_,*), intent(in)  plkcon,
integer, dimension(0:ntmat_,*), intent(in)  nplkcon,
real*8, dimension(27,mi(1),*), intent(in)  xstiff,
integer, intent(in)  npmat_,
real*8, intent(in)  dtime,
character*80, dimension(*), intent(in)  matname,
integer, dimension(*), intent(in)  mi,
integer, intent(in)  ncmat_,
integer, dimension(2), intent(in)  mass,
integer, intent(in)  stiffness,
integer, intent(in)  buckling,
integer, intent(in)  rhsi,
integer, intent(in)  intscheme,
real*8, dimension(*), intent(in)  physcon,
real*8, dimension(0:3,ntmat_,*), intent(in)  shcon,
integer, dimension(*), intent(in)  nshcon,
real*8, dimension(0:6,ntmat_,*), intent(in)  cocon,
integer, dimension(2,*), intent(in)  ncocon,
real*8, intent(in)  ttime,
real*8, intent(in)  time,
integer, intent(in)  istep,
integer, intent(in)  iinc,
integer, intent(in)  coriolis,
integer, dimension(3,*), intent(in)  ibody,
real*8, dimension(2,*), intent(in)  xloadold,
real*8, intent(in)  reltime,
real*8, dimension(0:mi(2),*), intent(in)  veold,
real*8, dimension(2,*), intent(inout)  springarea,
integer, intent(in)  nstate_,
real*8, dimension(nstate_,mi(1),*), intent(in)  xstateini,
real*8, dimension(nstate_,mi(1),*), intent(inout)  xstate,
real*8, dimension(mi(3),*), intent(in)  thicke,
integer, dimension(*), intent(in)  integerglob,
real*8, dimension(*), intent(in)  doubleglob,
character*81, dimension(3,*), intent(in)  tieset,
integer, dimension(*), intent(in)  istartset,
integer, dimension(*), intent(in)  iendset,
integer, dimension(*), intent(in)  ialset,
integer, intent(in)  ntie,
integer, intent(in)  nasym,
real*8, dimension(3,*), intent(in)  pslavsurf,
real*8, dimension(6,*), intent(in)  pmastsurf,
integer, intent(in)  mortar,
real*8, dimension(3,9,*), intent(in)  clearini,
integer, dimension(*), intent(in)  ielprop,
real*8, dimension(*), intent(in)  prop,
integer, intent(in)  ne0,
real*8, dimension(0:mi(2),*)  fnext,
integer, intent(in)  nea,
integer, intent(in)  neb,
integer  kscale,
integer, dimension(*)  iponoel,
integer, dimension(2,*)  inoel,
integer  network 
)
36 !
37 ! filling the stiffness matrix in spare matrix format (sm)
38 !
39  implicit none
40 !
41  integer mass(2),stiffness,buckling,rhsi,stiffonly(2),coriolis
42 !
43  character*8 lakon(*)
44  character*20 sideload(*)
45  character*80 matname(*)
46  character*81 tieset(3,*)
47 !
48  integer kon(*),nodeboun(*),ndirboun(*),ipompc(*),nodempc(3,*),
49  & nodeforc(2,*),ndirforc(*),nelemload(2,*),icol(*),jq(*),ikmpc(*),
50  & ilmpc(*),ikboun(*),ilboun(*),mi(*),nstate_,ne0,nasym,
51  & nactdof(0:mi(2),*),irow(*),icolumn,ialset(*),ielprop(*),
52  & nelcon(2,*),nrhcon(*),nalcon(2,*),ielmat(mi(3),*),ntie,
53  & ielorien(mi(3),*),integerglob(*),istartset(*),iendset(*),
54  & ipkon(*),intscheme,ncocon(2,*),nshcon(*),ipobody(2,*),nbody,
55  & ibody(3,*),nk,ne,nboun,nmpc,nforc,nload,neq(2),nzl,nmethod,
56  & ithermal(2),iprestr,iperturb(*),nzs(3),i,j,k,l,m,idist,jj,
57  & ll,id,id1,id2,ist,ist1,ist2,index,jdof1,jdof2,idof1,idof2,
58  & mpc1,mpc2,index1,index2,jdof,node1,node2,kflag,icalccg,
59  & ntmat_,indexe,nope,norien,iexpl,i0,ncmat_,istep,iinc,
60  & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),npmat_,mortar,
61  & nea,neb,kscale,iponoel(*),inoel(2,*),network,ndof
62 !
63  real*8 co(3,*),xboun(*),coefmpc(*),xforc(*),xload(2,*),p1(3),
64  & p2(3),ad(*),au(*),bodyf(3),fext(*),xloadold(2,*),reltime,
65  & t0(*),t1(*),prestr(6,mi(1),*),vold(0:mi(2),*),s(60,60),
66  & ff(60),fnext(0:mi(2),*),
67  & sti(6,mi(1),*),sm(60,60),stx(6,mi(1),*),adb(*),aub(*),
68  & elcon(0:ncmat_,ntmat_,*),rhcon(0:1,ntmat_,*),springarea(2,*),
69  & alcon(0:6,ntmat_,*),physcon(*),cocon(0:6,ntmat_,*),prop(*),
70  & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
71  & shcon(0:3,ntmat_,*),alzero(*),orab(7,*),xbody(7,*),cgr(4,*),
72  & plicon(0:2*npmat_,ntmat_,*),plkcon(0:2*npmat_,ntmat_,*),
73  & xstiff(27,mi(1),*),veold(0:mi(2),*),om,valu2,value,dtime,ttime,
74  & time,thicke(mi(3),*),doubleglob(*),clearini(3,9,*),
75  & pslavsurf(3,*),pmastsurf(6,*)
76 !
77  intent(in) co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,
78  & xboun,nboun,
79  & ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,xforc,
80  & nforc,nelemload,sideload,nload,xbody,ipobody,nbody,
81  & nactdof,icol,jq,irow,neq,nzl,
82  & ikmpc,ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,
83  & nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,
84  & t0,t1,ithermal,prestr,
85  & iprestr,vold,iperturb,sti,nzs,stx,iexpl,plicon,
86  & nplicon,plkcon,nplkcon,xstiff,npmat_,dtime,
87  & matname,mi,ncmat_,mass,stiffness,buckling,rhsi,intscheme,
88  & physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc,
89  & coriolis,ibody,xloadold,reltime,veold,nstate_,
90  & xstateini,thicke,integerglob,doubleglob,
91  & tieset,istartset,iendset,ialset,ntie,nasym,pslavsurf,pmastsurf,
92  & mortar,clearini,ielprop,prop,ne0,nea,neb
93 !
94  intent(inout) fext,ad,au,adb,aub,xload,nmethod,cgr,springarea,
95  & xstate
96 !
97  kflag=2
98  i0=0
99  icalccg=0
100 c write(*,*) loc(kflag)
101 c write(*,*) loc(s)
102 c write(*,*) loc(sm)
103 c write(*,*) loc(ff)
104 c write(*,*) loc(index1)
105 !
106  if((stiffness.eq.1).and.(mass(1).eq.0).and.(buckling.eq.0)) then
107  stiffonly(1)=1
108  else
109  stiffonly(1)=0
110  endif
111  if((stiffness.eq.1).and.(mass(2).eq.0).and.(buckling.eq.0)) then
112  stiffonly(2)=1
113  else
114  stiffonly(2)=0
115  endif
116 !
117  if(rhsi.eq.1) then
118 !
119 ! distributed forces (body forces or thermal loads or
120 ! residual stresses or distributed face loads)
121 !
122  if((nbody.ne.0).or.(ithermal(1).ne.0).or.
123  & (iprestr.ne.0).or.(nload.ne.0)) then
124  idist=1
125  else
126  idist=0
127  endif
128 !
129  endif
130 !
131  if((ithermal(1).le.1).or.(ithermal(1).eq.3)) then
132 !
133 ! mechanical analysis: loop over all elements
134 !
135  do i=nea,neb
136 !
137  if(ipkon(i).lt.0) cycle
138  indexe=ipkon(i)
139 c Bernhardi start
140  if(lakon(i)(1:5).eq.'C3D8I') then
141  nope=11
142  ndof=3
143  elseif(lakon(i)(4:5).eq.'20') then
144 c Bernhardi end
145  nope=20
146  ndof=3
147  elseif(lakon(i)(4:4).eq.'8') then
148  nope=8
149  ndof=3
150  elseif(lakon(i)(4:5).eq.'10') then
151  nope=10
152  ndof=3
153  elseif(lakon(i)(4:4).eq.'4') then
154  nope=4
155  ndof=3
156  elseif(lakon(i)(4:5).eq.'15') then
157  nope=15
158  ndof=3
159  elseif(lakon(i)(4:4).eq.'6') then
160  nope=6
161  ndof=3
162  elseif((lakon(i)(1:2).eq.'ES').and.(lakon(i)(7:7).ne.'F')) then
163 !
164 ! spring and contact spring elements (NO dashpot elements
165 ! = ED... elements)
166 !
167  nope=ichar(lakon(i)(8:8))-47
168  ndof=3
169 !
170 ! local contact spring number
171 ! if friction is involved, the contact spring element
172 ! matrices are determined in mafillsmas.f
173 !
174  if(lakon(i)(7:7).eq.'C') then
175  if(nasym.eq.1) cycle
176  if(mortar.eq.1) nope=kon(indexe)
177  endif
178  elseif(lakon(i)(1:4).eq.'MASS') then
179  nope=1
180  ndof=3
181  elseif(lakon(i)(1:1).eq.'U') then
182 c if(lakon(i)(7:7).eq.' ') then
183 c nope=ichar(lakon(i)(8:8))-48
184 c else
185 c nope=10*(ichar(lakon(i)(7:7))-48)
186 c & +ichar(lakon(i)(8:8))-48
187 c endif
188 c ndof=ichar(lakon(i)(6:6))-48
189  ndof=ichar(lakon(i)(7:7))
190  nope=ichar(lakon(i)(8:8))
191  else
192  cycle
193  endif
194 !
195  om=0.d0
196 !
197  if((nbody.gt.0).and.(lakon(i)(1:1).ne.'E')) then
198 !
199 ! assigning centrifugal forces
200 !
201  bodyf(1)=0.
202  bodyf(2)=0.
203  bodyf(3)=0.
204 !
205  index=i
206  do
207  j=ipobody(1,index)
208  if(j.eq.0) exit
209  if(ibody(1,j).eq.1) then
210  om=xbody(1,j)
211  p1(1)=xbody(2,j)
212  p1(2)=xbody(3,j)
213  p1(3)=xbody(4,j)
214  p2(1)=xbody(5,j)
215  p2(2)=xbody(6,j)
216  p2(3)=xbody(7,j)
217 !
218 ! assigning gravity forces
219 !
220  elseif(ibody(1,j).eq.2) then
221  bodyf(1)=bodyf(1)+xbody(1,j)*xbody(2,j)
222  bodyf(2)=bodyf(2)+xbody(1,j)*xbody(3,j)
223  bodyf(3)=bodyf(3)+xbody(1,j)*xbody(4,j)
224 !
225 ! assigning newton gravity forces
226 !
227  elseif(ibody(1,j).eq.3) then
228  call newton(icalccg,ne,ipkon,lakon,kon,t0,co,rhcon,
229  & nrhcon,ntmat_,physcon,i,cgr,bodyf,ielmat,ithermal,
230  & vold,mi)
231  endif
232  index=ipobody(2,index)
233  if(index.eq.0) exit
234  enddo
235  endif
236 !
237  if(lakon(i)(1:1).ne.'U') then
238  call e_c3d(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,ff,i,
239  & nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
240  & alzero,ielmat,ielorien,norien,orab,ntmat_,
241  & t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload,
242  & nload,idist,sti,stx,iexpl,plicon,
243  & nplicon,plkcon,nplkcon,xstiff,npmat_,
244  & dtime,matname,mi(1),ncmat_,mass(1),stiffness,buckling,
245  & rhsi,intscheme,ttime,time,istep,iinc,coriolis,xloadold,
246  & reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
247  & springarea,nstate_,xstateini,xstate,ne0,ipkon,thicke,
248  & integerglob,doubleglob,tieset,istartset,
249  & iendset,ialset,ntie,nasym,pslavsurf,pmastsurf,mortar,
250  & clearini,ielprop,prop,kscale)
251  else
252  call e_c3d_u(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,ff,i,
253  & nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
254  & alzero,ielmat,ielorien,norien,orab,ntmat_,
255  & t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload,
256  & nload,idist,sti,stx,iexpl,plicon,
257  & nplicon,plkcon,nplkcon,xstiff,npmat_,
258  & dtime,matname,mi(1),ncmat_,mass(1),stiffness,buckling,
259  & rhsi,intscheme,ttime,time,istep,iinc,coriolis,xloadold,
260  & reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
261  & ne0,ipkon,thicke,
262  & integerglob,doubleglob,tieset,istartset,
263  & iendset,ialset,ntie,nasym,
264  & ielprop,prop)
265  endif
266 !
267  do jj=1,ndof*nope
268 !
269  j=(jj-1)/ndof+1
270  k=jj-ndof*(j-1)
271 !
272  node1=kon(indexe+j)
273  jdof1=nactdof(k,node1)
274 !
275  do ll=jj,ndof*nope
276 !
277  l=(ll-1)/ndof+1
278  m=ll-ndof*(l-1)
279 !
280  node2=kon(indexe+l)
281  jdof2=nactdof(m,node2)
282 !
283 ! check whether one of the DOF belongs to a SPC or MPC
284 !
285  if((jdof1.gt.0).and.(jdof2.gt.0)) then
286  if(stiffonly(1).eq.1) then
287  call add_sm_st(au,ad,jq,irow,jdof1,jdof2,
288  & s(jj,ll),jj,ll)
289  else
290  call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1,jdof2,
291  & s(jj,ll),sm(jj,ll),jj,ll)
292  endif
293  elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
294 !
295 ! idof1: genuine DOF
296 ! idof2: nominal DOF of the SPC/MPC
297 !
298  if(jdof1.le.0) then
299  idof1=jdof2
300  idof2=jdof1
301  else
302  idof1=jdof1
303  idof2=jdof2
304  endif
305  if(nmpc.gt.0) then
306  if(idof2.ne.2*(idof2/2)) then
307 !
308 ! regular DOF / MPC
309 !
310  id=(-idof2+1)/2
311  ist=ipompc(id)
312  index=nodempc(3,ist)
313  if(index.eq.0) cycle
314  do
315  idof2=nactdof(nodempc(2,index),nodempc(1,index))
316  value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
317  if(idof1.eq.idof2) value=2.d0*value
318  if(idof2.gt.0) then
319  if(stiffonly(1).eq.1) then
320  call add_sm_st(au,ad,jq,irow,idof1,
321  & idof2,value,i0,i0)
322  else
323  valu2=-coefmpc(index)*sm(jj,ll)/
324  & coefmpc(ist)
325 !
326  if(idof1.eq.idof2) valu2=2.d0*valu2
327 !
328  call add_sm_ei(au,ad,aub,adb,jq,irow,
329  & idof1,idof2,value,valu2,i0,i0)
330  endif
331  endif
332  index=nodempc(3,index)
333  if(index.eq.0) exit
334  enddo
335  cycle
336  endif
337  endif
338 !
339 ! regular DOF / SPC
340 !
341  if(rhsi.eq.1) then
342  elseif(nmethod.eq.2) then
343  value=s(jj,ll)
344  icolumn=neq(2)-idof2/2
345  call add_bo_st(au,jq,irow,idof1,icolumn,value)
346  endif
347  else
348  idof1=jdof1
349  idof2=jdof2
350  mpc1=0
351  mpc2=0
352  if(nmpc.gt.0) then
353  if(idof1.ne.2*(idof1/2)) mpc1=1
354  if(idof2.ne.2*(idof2/2)) mpc2=1
355  endif
356  if((mpc1.eq.1).and.(mpc2.eq.1)) then
357  id1=(-idof1+1)/2
358  id2=(-idof2+1)/2
359  if(id1.eq.id2) then
360 !
361 ! MPC id1 / MPC id1
362 !
363  ist=ipompc(id1)
364  index1=nodempc(3,ist)
365  if(index1.eq.0) cycle
366  do
367  idof1=nactdof(nodempc(2,index1),
368  & nodempc(1,index1))
369  index2=index1
370  do
371  idof2=nactdof(nodempc(2,index2),
372  & nodempc(1,index2))
373  value=coefmpc(index1)*coefmpc(index2)*
374  & s(jj,ll)/coefmpc(ist)/coefmpc(ist)
375  if((idof1.gt.0).and.(idof2.gt.0)) then
376  if(stiffonly(1).eq.1) then
377  call add_sm_st(au,ad,jq,irow,
378  & idof1,idof2,value,i0,i0)
379  else
380  valu2=coefmpc(index1)*coefmpc(index2)*
381  & sm(jj,ll)/coefmpc(ist)/coefmpc(ist)
382  call add_sm_ei(au,ad,aub,adb,jq,
383  & irow,idof1,idof2,value,valu2,i0,i0)
384  endif
385  endif
386 !
387  index2=nodempc(3,index2)
388  if(index2.eq.0) exit
389  enddo
390  index1=nodempc(3,index1)
391  if(index1.eq.0) exit
392  enddo
393  else
394 !
395 ! MPC id1 / MPC id2
396 !
397  ist1=ipompc(id1)
398  index1=nodempc(3,ist1)
399  if(index1.eq.0) cycle
400  do
401  idof1=nactdof(nodempc(2,index1),
402  & nodempc(1,index1))
403  ist2=ipompc(id2)
404  index2=nodempc(3,ist2)
405  if(index2.eq.0) then
406  index1=nodempc(3,index1)
407  if(index1.eq.0) then
408  exit
409  else
410  cycle
411  endif
412  endif
413  do
414  idof2=nactdof(nodempc(2,index2),
415  & nodempc(1,index2))
416  value=coefmpc(index1)*coefmpc(index2)*
417  & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
418  if(idof1.eq.idof2) value=2.d0*value
419  if((idof1.gt.0).and.(idof2.gt.0)) then
420  if(stiffonly(1).eq.1) then
421  call add_sm_st(au,ad,jq,irow,
422  & idof1,idof2,value,i0,i0)
423  else
424  valu2=coefmpc(index1)*coefmpc(index2)*
425  & sm(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
426 !
427  if(idof1.eq.idof2) valu2=2.d0*valu2
428 !
429  call add_sm_ei(au,ad,aub,adb,jq,
430  & irow,idof1,idof2,value,valu2,i0,i0)
431  endif
432  endif
433 !
434  index2=nodempc(3,index2)
435  if(index2.eq.0) exit
436  enddo
437  index1=nodempc(3,index1)
438  if(index1.eq.0) exit
439  enddo
440  endif
441  endif
442  endif
443  enddo
444 !
445  if(rhsi.eq.1) then
446 !
447 ! distributed forces
448 !
449  if(idist.ne.0) then
450 !
451 ! updating the external force vector for dynamic
452 ! calculations
453 !
454  if(nmethod.eq.4) fnext(k,node1)=fnext(k,node1)+ff(jj)
455 !
456  if(jdof1.le.0) then
457  if(nmpc.ne.0) then
458  idof1=jdof1
459  if(idof1.ne.2*(idof1/2)) then
460  id=(-idof1+1)/2
461  ist=ipompc(id)
462  index=nodempc(3,ist)
463  if(index.eq.0) cycle
464  do
465  jdof1=nactdof(nodempc(2,index),
466  & nodempc(1,index))
467  if(jdof1.gt.0) then
468  fext(jdof1)=fext(jdof1)
469  & -coefmpc(index)*ff(jj)
470  & /coefmpc(ist)
471  endif
472  index=nodempc(3,index)
473  if(index.eq.0) exit
474  enddo
475  endif
476  endif
477  cycle
478  endif
479  fext(jdof1)=fext(jdof1)+ff(jj)
480  endif
481  endif
482 !
483  enddo
484  enddo
485 !
486  endif
487  if(ithermal(1).gt.1) then
488 !
489 ! thermal analysis: loop over all elements
490 !
491  do i=nea,neb
492 !
493  if(ipkon(i).lt.0) cycle
494  indexe=ipkon(i)
495  if(lakon(i)(4:5).eq.'20') then
496  nope=20
497  elseif(lakon(i)(4:4).eq.'8') then
498  nope=8
499  elseif(lakon(i)(4:5).eq.'10') then
500  nope=10
501  elseif(lakon(i)(4:4).eq.'4') then
502  nope=4
503  elseif(lakon(i)(4:5).eq.'15') then
504  nope=15
505  elseif(lakon(i)(4:4).eq.'6') then
506  nope=6
507  elseif((lakon(i)(1:1).eq.'E').and.(lakon(i)(7:7).ne.'A')) then
508 !
509 ! contact spring and advection elements
510 !
511  nope=ichar(lakon(i)(8:8))-47
512 !
513 ! local contact spring number
514 !
515  if(lakon(i)(7:7).eq.'C') then
516  if(mortar.eq.1) nope=kon(indexe)
517  endif
518  elseif((lakon(i)(1:2).eq.'D ').or.
519  & ((lakon(i)(1:1).eq.'D').and.(network.eq.1))) then
520 !
521 ! asymmetrical contribution -> mafillsmas.f
522 !
523  cycle
524  else
525  cycle
526  endif
527 !
528  call e_c3d_th(co,nk,kon,lakon(i),s,sm,
529  & ff,i,nmethod,rhcon,nrhcon,ielmat,ielorien,norien,orab,
530  & ntmat_,t0,t1,ithermal,vold,iperturb,nelemload,
531  & sideload,xload,nload,idist,iexpl,dtime,
532  & matname,mi(1),mass(2),stiffness,buckling,rhsi,intscheme,
533  & physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc,
534  & xstiff,xloadold,reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,
535  & ilmpc,springarea,plkcon,nplkcon,npmat_,ncmat_,elcon,nelcon,
536  & lakon,pslavsurf,pmastsurf,mortar,clearini,plicon,nplicon,
537  & ipkon,ielprop,prop,iponoel,inoel,sti,xstateini,xstate,
538  & nstate_,network,ipobody,xbody,ibody)
539 !
540  do jj=1,nope
541 !
542  j=jj
543 !
544  node1=kon(indexe+j)
545  jdof1=nactdof(0,node1)
546 !
547  do ll=jj,nope
548 !
549  l=ll
550 !
551  node2=kon(indexe+l)
552  jdof2=nactdof(0,node2)
553 !
554 ! check whether one of the DOF belongs to a SPC or MPC
555 !
556  if((jdof1.gt.0).and.(jdof2.gt.0)) then
557  if(stiffonly(2).eq.1) then
558  call add_sm_st(au,ad,jq,irow,jdof1,jdof2,
559  & s(jj,ll),jj,ll)
560  else
561  call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1,jdof2,
562  & s(jj,ll),sm(jj,ll),jj,ll)
563  endif
564  elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
565 !
566 ! idof1: genuine DOF
567 ! idof2: nominal DOF of the SPC/MPC
568 !
569  if(jdof1.le.0) then
570  idof1=jdof2
571  idof2=jdof1
572  else
573  idof1=jdof1
574  idof2=jdof2
575  endif
576  if(nmpc.gt.0) then
577  if(idof2.ne.2*(idof2/2)) then
578 !
579 ! regular DOF / MPC
580 !
581  id=(-idof2+1)/2
582  ist=ipompc(id)
583  index=nodempc(3,ist)
584  if(index.eq.0) cycle
585  do
586  idof2=nactdof(nodempc(2,index),nodempc(1,index))
587  value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
588  if(idof1.eq.idof2) value=2.d0*value
589  if(idof2.gt.0) then
590  if(stiffonly(2).eq.1) then
591  call add_sm_st(au,ad,jq,irow,idof1,
592  & idof2,value,i0,i0)
593  else
594  valu2=-coefmpc(index)*sm(jj,ll)/
595  & coefmpc(ist)
596 !
597  if(idof1.eq.idof2) valu2=2.d0*valu2
598 !
599  call add_sm_ei(au,ad,aub,adb,jq,irow,
600  & idof1,idof2,value,valu2,i0,i0)
601  endif
602  endif
603  index=nodempc(3,index)
604  if(index.eq.0) exit
605  enddo
606  cycle
607  endif
608  endif
609 !
610 ! regular DOF / SPC
611 !
612  if(rhsi.eq.1) then
613  elseif(nmethod.eq.2) then
614  value=s(jj,ll)
615  icolumn=neq(2)-idof2/2
616  call add_bo_st(au,jq,irow,idof1,icolumn,value)
617  endif
618  else
619  idof1=jdof1
620  idof2=jdof2
621  mpc1=0
622  mpc2=0
623  if(nmpc.gt.0) then
624  if(idof1.ne.2*(idof1/2)) mpc1=1
625  if(idof2.ne.2*(idof2/2)) mpc2=1
626  endif
627  if((mpc1.eq.1).and.(mpc2.eq.1)) then
628  id1=(-idof1+1)/2
629  id2=(-idof2+1)/2
630  if(id1.eq.id2) then
631 !
632 ! MPC id1 / MPC id1
633 !
634  ist=ipompc(id1)
635  index1=nodempc(3,ist)
636  if(index1.eq.0) cycle
637  do
638  idof1=nactdof(nodempc(2,index1),
639  & nodempc(1,index1))
640  index2=index1
641  do
642  idof2=nactdof(nodempc(2,index2),
643  & nodempc(1,index2))
644  value=coefmpc(index1)*coefmpc(index2)*
645  & s(jj,ll)/coefmpc(ist)/coefmpc(ist)
646  if((idof1.gt.0).and.(idof2.gt.0)) then
647  if(stiffonly(2).eq.1) then
648  call add_sm_st(au,ad,jq,irow,
649  & idof1,idof2,value,i0,i0)
650  else
651  valu2=coefmpc(index1)*coefmpc(index2)*
652  & sm(jj,ll)/coefmpc(ist)/coefmpc(ist)
653  call add_sm_ei(au,ad,aub,adb,jq,
654  & irow,idof1,idof2,value,valu2,i0,i0)
655  endif
656  endif
657 !
658  index2=nodempc(3,index2)
659  if(index2.eq.0) exit
660  enddo
661  index1=nodempc(3,index1)
662  if(index1.eq.0) exit
663  enddo
664  else
665 !
666 ! MPC id1 / MPC id2
667 !
668  ist1=ipompc(id1)
669  index1=nodempc(3,ist1)
670  if(index1.eq.0) cycle
671  do
672  idof1=nactdof(nodempc(2,index1),
673  & nodempc(1,index1))
674  ist2=ipompc(id2)
675  index2=nodempc(3,ist2)
676  if(index2.eq.0) then
677  index1=nodempc(3,index1)
678  if(index1.eq.0) then
679  exit
680  else
681  cycle
682  endif
683  endif
684  do
685  idof2=nactdof(nodempc(2,index2),
686  & nodempc(1,index2))
687  value=coefmpc(index1)*coefmpc(index2)*
688  & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
689  if(idof1.eq.idof2) value=2.d0*value
690  if((idof1.gt.0).and.(idof2.gt.0)) then
691  if(stiffonly(2).eq.1) then
692  call add_sm_st(au,ad,jq,irow,
693  & idof1,idof2,value,i0,i0)
694  else
695  valu2=coefmpc(index1)*coefmpc(index2)*
696  & sm(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
697 !
698  if(idof1.eq.idof2) valu2=2.d0*valu2
699 !
700  call add_sm_ei(au,ad,aub,adb,jq,
701  & irow,idof1,idof2,value,valu2,i0,i0)
702  endif
703  endif
704 !
705  index2=nodempc(3,index2)
706  if(index2.eq.0) exit
707  enddo
708  index1=nodempc(3,index1)
709  if(index1.eq.0) exit
710  enddo
711  endif
712  endif
713  endif
714  enddo
715 !
716  if(rhsi.eq.1) then
717 !
718 ! distributed forces
719 !
720  if(idist.ne.0) then
721  if(jdof1.le.0) then
722  if(nmpc.ne.0) then
723  idof1=jdof1
724  if(idof1.ne.2*(idof1/2)) then
725  id=(-idof1+1)/2
726  ist=ipompc(id)
727  index=nodempc(3,ist)
728  if(index.eq.0) cycle
729  do
730  jdof1=nactdof(nodempc(2,index),
731  & nodempc(1,index))
732  if(jdof1.gt.0) then
733  fext(jdof1)=fext(jdof1)
734  & -coefmpc(index)*ff(jj)
735  & /coefmpc(ist)
736  endif
737  index=nodempc(3,index)
738  if(index.eq.0) exit
739  enddo
740  endif
741  endif
742  cycle
743  endif
744  fext(jdof1)=fext(jdof1)+ff(jj)
745  endif
746  endif
747 !
748  enddo
749  enddo
750 !
751  endif
752 !
753  return
newton
subroutine newton(icalccg, ne, ipkon, lakon, kon, t0, co, rhcon, nrhcon, ntmat_, physcon, nelem, cgr, bodyf, ielmat, ithermal, vold, mi)
Definition: newton.f:22
e_c3d
subroutine e_c3d(co, kon, lakonl, p1, p2, omx, bodyfx, nbody, s, sm, ff, nelem, nmethod, 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_, mass, stiffness, buckling, rhsi, intscheme, ttime, time, istep, iinc, coriolis, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, veold, springarea, nstate_, xstateini, xstate, ne0, ipkon, thicke, integerglob, doubleglob, tieset, istartset, iendset, ialset, ntie, nasym, pslavsurf, pmastsurf, mortar, clearini, ielprop, prop, kscale)
Definition: e_c3d.f:31
add_bo_st
subroutine add_bo_st(au, jq, irow, i, j, value)
Definition: add_bo_st.f:20
e_c3d_u
subroutine e_c3d_u(co, kon, lakonl, p1, p2, omx, bodyfx, nbody, s, sm, ff, nelem, nmethod, 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_, mass, stiffness, buckling, rhsi, intscheme, ttime, time, istep, iinc, coriolis, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, veold, ne0, ipkon, thicke, integerglob, doubleglob, tieset, istartset, iendset, ialset, ntie, nasym, ielprop, prop)
Definition: e_c3d_u.f:31
add_sm_ei
subroutine add_sm_ei(au, ad, aub, adb, jq, irow, i, j, value, valuem, i0, i1)
Definition: add_sm_ei.f:21
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
e_c3d_th
subroutine e_c3d_th(co, nk, kon, lakonl, s, sm, ff, nelem, nmethod, rhcon, nrhcon, ielmat, ielorien, norien, orab, ntmat_, t0, t1, ithermal, vold, iperturb, nelemload, sideload, xload, nload, idist, iexpl, dtime, matname, mi, mass, stiffness, buckling, rhsi, intscheme, physcon, shcon, nshcon, cocon, ncocon, ttime, time, istep, iinc, xstiff, xloadold, reltime, ipompc, nodempc, coefmpc, nmpc, ikmpc, ilmpc, springarea, plkcon, nplkcon, npmat_, ncmat_, elcon, nelcon, lakon, pslavsurf, pmastsurf, mortar, clearini, plicon, nplicon, ipkon, ielprop, prop, iponoel, inoel, sti, xstateini, xstate, nstate_, network, ipobody, xbody, ibody)
Definition: e_c3d_th.f:31
Hosted by OpenAircraft.com, (Michigan UAV, LLC)