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

Go to the source code of this file.

Functions/Subroutines

subroutine mafillsm_company (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, kscale, iponoel, inoel, network)
 

Function/Subroutine Documentation

◆ mafillsm_company()

subroutine mafillsm_company ( 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(inout)  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  kscale,
integer, dimension(*), intent(in)  iponoel,
integer, dimension(2,*), intent(in)  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  & kscale,iponoel(*),inoel(2,*),network
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,
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,iponoel,inoel
93 !
94  intent(inout) fext,ad,au,adb,aub,xload,nmethod,cgr,springarea,
95  & xstate,nzl
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 ! determining nzl
118 !
119  nzl=0
120  do i=neq(2),1,-1
121  if(icol(i).gt.0) then
122  nzl=i
123  exit
124  endif
125  enddo
126 !
127 ! initializing the matrices
128 !
129  if(buckling.ne.1) then
130  do i=1,neq(2)
131  ad(i)=0.d0
132  enddo
133  do i=1,nzs(3)
134  au(i)=0.d0
135  enddo
136  endif
137 !
138  if(rhsi.eq.1) then
139  do i=1,neq(2)
140  fext(i)=0.d0
141  enddo
142  endif
143  if((mass(1).eq.1).or.(buckling.eq.1)) then
144  do i=1,neq(1)
145  adb(i)=0.d0
146  enddo
147  do i=1,nzs(1)
148  aub(i)=0.d0
149  enddo
150  endif
151  if(mass(2).eq.1) then
152  do i=neq(1)+1,neq(2)
153  adb(i)=0.d0
154  enddo
155  do i=nzs(1)+1,nzs(2)
156  aub(i)=0.d0
157  enddo
158  endif
159 !
160  if(rhsi.eq.1) then
161 !
162 ! distributed forces (body forces or thermal loads or
163 ! residual stresses or distributed face loads)
164 !
165  if((nbody.ne.0).or.(ithermal(1).ne.0).or.
166  & (iprestr.ne.0).or.(nload.ne.0)) then
167  idist=1
168  else
169  idist=0
170  endif
171 !
172  endif
173 !
174  if((ithermal(1).le.1).or.(ithermal(1).eq.3)) then
175 !
176 ! mechanical analysis: loop over all elements
177 !
178  do i=1,ne
179 !
180  if(ipkon(i).lt.0) cycle
181  indexe=ipkon(i)
182 c Bernhardi start
183  if(lakon(i)(1:5).eq.'C3D8I') then
184  nope=11
185  elseif(lakon(i)(4:5).eq.'20') then
186 c Bernhardi end
187  nope=20
188  elseif(lakon(i)(4:4).eq.'8') then
189  nope=8
190  elseif(lakon(i)(4:5).eq.'10') then
191  nope=10
192  elseif(lakon(i)(4:4).eq.'4') then
193  nope=4
194  elseif(lakon(i)(4:5).eq.'15') then
195  nope=15
196  elseif(lakon(i)(4:4).eq.'6') then
197  nope=6
198  elseif((lakon(i)(1:2).eq.'ES').and.(lakon(i)(7:7).ne.'F')) then
199 !
200 ! spring and contact spring elements (NO dashpot elements
201 ! = ED... elements)
202 !
203 c read(lakon(i)(8:8),'(i1)') nope
204  nope=ichar(lakon(i)(8:8))-47
205 c nope=nope+1
206 !
207 ! local contact spring number
208 ! if friction is involved, the contact spring element
209 ! matrices are determined in mafillsmas.f
210 !
211  if(lakon(i)(7:7).eq.'C') then
212  if(nasym.eq.1) cycle
213  if(mortar.eq.1) nope=kon(indexe)
214  endif
215  else
216  cycle
217  endif
218 !
219  om=0.d0
220 !
221  if((nbody.gt.0).and.(lakon(i)(1:1).ne.'E')) then
222 !
223 ! assigning centrifugal forces
224 !
225  bodyf(1)=0.
226  bodyf(2)=0.
227  bodyf(3)=0.
228 !
229  index=i
230  do
231  j=ipobody(1,index)
232  if(j.eq.0) exit
233  if(ibody(1,j).eq.1) then
234  om=xbody(1,j)
235  p1(1)=xbody(2,j)
236  p1(2)=xbody(3,j)
237  p1(3)=xbody(4,j)
238  p2(1)=xbody(5,j)
239  p2(2)=xbody(6,j)
240  p2(3)=xbody(7,j)
241 !
242 ! assigning gravity forces
243 !
244  elseif(ibody(1,j).eq.2) then
245  bodyf(1)=bodyf(1)+xbody(1,j)*xbody(2,j)
246  bodyf(2)=bodyf(2)+xbody(1,j)*xbody(3,j)
247  bodyf(3)=bodyf(3)+xbody(1,j)*xbody(4,j)
248 !
249 ! assigning newton gravity forces
250 !
251  elseif(ibody(1,j).eq.3) then
252  call newton(icalccg,ne,ipkon,lakon,kon,t0,co,rhcon,
253  & nrhcon,ntmat_,physcon,i,cgr,bodyf,ielmat,ithermal,
254  & vold,mi)
255  endif
256  index=ipobody(2,index)
257  if(index.eq.0) exit
258  enddo
259  endif
260 !
261  call e_c3d(co,kon,lakon(i),p1,p2,om,bodyf,nbody,s,sm,ff,i,
262  & nmethod,elcon,nelcon,rhcon,nrhcon,alcon,nalcon,
263  & alzero,ielmat,ielorien,norien,orab,ntmat_,
264  & t0,t1,ithermal,vold,iperturb,nelemload,sideload,xload,
265  & nload,idist,sti,stx,iexpl,plicon,
266  & nplicon,plkcon,nplkcon,xstiff,npmat_,
267  & dtime,matname,mi(1),ncmat_,mass(1),stiffness,buckling,
268  & rhsi,intscheme,ttime,time,istep,iinc,coriolis,xloadold,
269  & reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,veold,
270  & springarea,nstate_,xstateini,xstate,ne0,ipkon,thicke,
271  & integerglob,doubleglob,tieset,istartset,
272  & iendset,ialset,ntie,nasym,pslavsurf,pmastsurf,mortar,
273  & clearini,ielprop,prop,kscale)
274 !
275  do jj=1,3*nope
276 !
277  j=(jj-1)/3+1
278  k=jj-3*(j-1)
279 !
280  node1=kon(indexe+j)
281  jdof1=nactdof(k,node1)
282 !
283  do ll=jj,3*nope
284 !
285  l=(ll-1)/3+1
286  m=ll-3*(l-1)
287 !
288  node2=kon(indexe+l)
289  jdof2=nactdof(m,node2)
290 !
291 ! check whether one of the DOF belongs to a SPC or MPC
292 !
293  if((jdof1.gt.0).and.(jdof2.gt.0)) then
294  if(stiffonly(1).eq.1) then
295  call add_sm_st(au,ad,jq,irow,jdof1,jdof2,
296  & s(jj,ll),jj,ll)
297  else
298  call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1,jdof2,
299  & s(jj,ll),sm(jj,ll),jj,ll)
300  endif
301  elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
302 !
303 ! idof1: genuine DOF
304 ! idof2: nominal DOF of the SPC/MPC
305 !
306  if(jdof1.le.0) then
307  idof1=jdof2
308 c idof2=(node1-1)*8+k
309  idof2=jdof1
310  else
311  idof1=jdof1
312 c idof2=(node2-1)*8+m
313  idof2=jdof2
314  endif
315  if(nmpc.gt.0) then
316 c call nident(ikmpc,idof2,nmpc,id)
317 c if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
318  if(idof2.ne.2*(idof2/2)) then
319 !
320 ! regular DOF / MPC
321 !
322 c id=ilmpc(id)
323  id=(-idof2+1)/2
324  ist=ipompc(id)
325  index=nodempc(3,ist)
326  if(index.eq.0) cycle
327  do
328  idof2=nactdof(nodempc(2,index),nodempc(1,index))
329  value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
330  if(idof1.eq.idof2) value=2.d0*value
331  if(idof2.gt.0) then
332  if(stiffonly(1).eq.1) then
333  call add_sm_st(au,ad,jq,irow,idof1,
334  & idof2,value,i0,i0)
335  else
336  valu2=-coefmpc(index)*sm(jj,ll)/
337  & coefmpc(ist)
338 c
339  if(idof1.eq.idof2) valu2=2.d0*valu2
340 c
341  call add_sm_ei(au,ad,aub,adb,jq,irow,
342  & idof1,idof2,value,valu2,i0,i0)
343  endif
344  endif
345  index=nodempc(3,index)
346  if(index.eq.0) exit
347  enddo
348  cycle
349  endif
350  endif
351 !
352 ! regular DOF / SPC
353 !
354  if(rhsi.eq.1) then
355  elseif(nmethod.eq.2) then
356  value=s(jj,ll)
357 c call nident(ikboun,idof2,nboun,id)
358 c icolumn=neq(2)+ilboun(id)
359  icolumn=neq(2)-idof2/2
360  call add_bo_st(au,jq,irow,idof1,icolumn,value)
361  endif
362  else
363 c idof1=(node1-1)*8+k
364 c idof2=(node2-1)*8+m
365  idof1=jdof1
366  idof2=jdof2
367  mpc1=0
368  mpc2=0
369  if(nmpc.gt.0) then
370 c call nident(ikmpc,idof1,nmpc,id1)
371 c if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
372 c call nident(ikmpc,idof2,nmpc,id2)
373 c if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
374  if(idof1.ne.2*(idof1/2)) mpc1=1
375  if(idof2.ne.2*(idof2/2)) mpc2=1
376  endif
377  if((mpc1.eq.1).and.(mpc2.eq.1)) then
378 c id1=ilmpc(id1)
379 c id2=ilmpc(id2)
380  id1=(-idof1+1)/2
381  id2=(-idof2+1)/2
382  if(id1.eq.id2) then
383 !
384 ! MPC id1 / MPC id1
385 !
386  ist=ipompc(id1)
387  index1=nodempc(3,ist)
388  if(index1.eq.0) cycle
389  do
390  idof1=nactdof(nodempc(2,index1),
391  & nodempc(1,index1))
392  index2=index1
393  do
394  idof2=nactdof(nodempc(2,index2),
395  & nodempc(1,index2))
396  value=coefmpc(index1)*coefmpc(index2)*
397  & s(jj,ll)/coefmpc(ist)/coefmpc(ist)
398  if((idof1.gt.0).and.(idof2.gt.0)) then
399  if(stiffonly(1).eq.1) then
400  call add_sm_st(au,ad,jq,irow,
401  & idof1,idof2,value,i0,i0)
402  else
403  valu2=coefmpc(index1)*coefmpc(index2)*
404  & sm(jj,ll)/coefmpc(ist)/coefmpc(ist)
405  call add_sm_ei(au,ad,aub,adb,jq,
406  & irow,idof1,idof2,value,valu2,i0,i0)
407  endif
408  endif
409 !
410  index2=nodempc(3,index2)
411  if(index2.eq.0) exit
412  enddo
413  index1=nodempc(3,index1)
414  if(index1.eq.0) exit
415  enddo
416  else
417 !
418 ! MPC id1 / MPC id2
419 !
420  ist1=ipompc(id1)
421  index1=nodempc(3,ist1)
422  if(index1.eq.0) cycle
423  do
424  idof1=nactdof(nodempc(2,index1),
425  & nodempc(1,index1))
426  ist2=ipompc(id2)
427  index2=nodempc(3,ist2)
428  if(index2.eq.0) then
429  index1=nodempc(3,index1)
430  if(index1.eq.0) then
431  exit
432  else
433  cycle
434  endif
435  endif
436  do
437  idof2=nactdof(nodempc(2,index2),
438  & nodempc(1,index2))
439  value=coefmpc(index1)*coefmpc(index2)*
440  & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
441  if(idof1.eq.idof2) value=2.d0*value
442  if((idof1.gt.0).and.(idof2.gt.0)) then
443  if(stiffonly(1).eq.1) then
444  call add_sm_st(au,ad,jq,irow,
445  & idof1,idof2,value,i0,i0)
446  else
447  valu2=coefmpc(index1)*coefmpc(index2)*
448  & sm(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
449 c
450  if(idof1.eq.idof2) valu2=2.d0*valu2
451 c
452  call add_sm_ei(au,ad,aub,adb,jq,
453  & irow,idof1,idof2,value,valu2,i0,i0)
454  endif
455  endif
456 !
457  index2=nodempc(3,index2)
458  if(index2.eq.0) exit
459  enddo
460  index1=nodempc(3,index1)
461  if(index1.eq.0) exit
462  enddo
463  endif
464 c elseif(((mpc1.eq.1).or.(mpc2.eq.1)).and.rhsi)
465 c & then
466 c if(mpc1.eq.1) then
467 c!
468 c! MPC id1 / SPC
469 c!
470 c call nident(ikboun,idof2,nboun,id2)
471 c idof2=ilboun(id2)
472 c ist1=ipompc(id1)
473 c index1=nodempc(3,ist1)
474 c if(index1.eq.0) cycle
475 c elseif(mpc2.eq.1) then
476 c!
477 c! MPC id2 / SPC
478 c!
479 c call nident(ikboun,idof1,nboun,id1)
480 c idof1=ilboun(id1)
481 c ist2=ipompc(id2)
482 c index2=nodempc(3,ist2)
483 c if(index2.eq.0) cycle
484 c endif
485  endif
486  endif
487  enddo
488 !
489  if(rhsi.eq.1) then
490 !
491 ! distributed forces
492 !
493  if(idist.ne.0) then
494 !
495 ! updating the external force vector for dynamic
496 ! calculations
497 !
498  if(nmethod.eq.4) fnext(k,node1)=fnext(k,node1)+ff(jj)
499 !
500  if(jdof1.le.0) then
501  if(nmpc.ne.0) then
502 c idof1=(node1-1)*8+k
503  idof1=jdof1
504 c call nident(ikmpc,idof1,nmpc,id)
505 c if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
506  if(idof1.ne.2*(idof1/2)) then
507 c id=ilmpc(id)
508  id=(-idof1+1)/2
509  ist=ipompc(id)
510  index=nodempc(3,ist)
511  if(index.eq.0) cycle
512  do
513  jdof1=nactdof(nodempc(2,index),
514  & nodempc(1,index))
515  if(jdof1.gt.0) then
516  fext(jdof1)=fext(jdof1)
517  & -coefmpc(index)*ff(jj)
518  & /coefmpc(ist)
519  endif
520  index=nodempc(3,index)
521  if(index.eq.0) exit
522  enddo
523  endif
524  endif
525  cycle
526  endif
527  fext(jdof1)=fext(jdof1)+ff(jj)
528  endif
529  endif
530 !
531  enddo
532  enddo
533 !
534  endif
535  if(ithermal(1).gt.1) then
536 !
537 ! thermal analysis: loop over all elements
538 !
539  do i=1,ne
540 !
541  if(ipkon(i).lt.0) cycle
542  indexe=ipkon(i)
543  if(lakon(i)(4:5).eq.'20') then
544  nope=20
545  elseif(lakon(i)(4:4).eq.'8') then
546  nope=8
547  elseif(lakon(i)(4:5).eq.'10') then
548  nope=10
549  elseif(lakon(i)(4:4).eq.'4') then
550  nope=4
551  elseif(lakon(i)(4:5).eq.'15') then
552  nope=15
553  elseif(lakon(i)(4:4).eq.'6') then
554  nope=6
555  elseif((lakon(i)(1:1).eq.'E').and.(lakon(i)(7:7).ne.'A')) then
556 !
557 ! contact spring and advection elements
558 !
559 c read(lakon(i)(8:8),'(i1)') nope
560  nope=ichar(lakon(i)(8:8))-47
561 c nope=nope+1
562 !
563 ! local contact spring number
564 !
565  if(lakon(i)(7:7).eq.'C') then
566  if(mortar.eq.1) nope=kon(indexe)
567  endif
568  elseif((lakon(i)(1:2).eq.'D ').or.
569  & ((lakon(i)(1:1).eq.'D').and.(network.eq.1))) then
570 !
571 ! asymmetrical contribution -> mafillsmas.f
572 !
573  cycle
574  else
575  cycle
576  endif
577 !
578  call e_c3d_th(co,nk,kon,lakon(i),s,sm,
579  & ff,i,nmethod,rhcon,nrhcon,ielmat,ielorien,norien,orab,
580  & ntmat_,t0,t1,ithermal,vold,iperturb,nelemload,
581  & sideload,xload,nload,idist,iexpl,dtime,
582  & matname,mi(1),mass(2),stiffness,buckling,rhsi,intscheme,
583  & physcon,shcon,nshcon,cocon,ncocon,ttime,time,istep,iinc,
584  & xstiff,xloadold,reltime,ipompc,nodempc,coefmpc,nmpc,ikmpc,
585  & ilmpc,springarea,plkcon,nplkcon,npmat_,ncmat_,elcon,nelcon,
586  & lakon,pslavsurf,pmastsurf,mortar,clearini,plicon,nplicon,
587  & ipkon,ielprop,prop,iponoel,inoel,sti,xstateini,xstate,
588  & nstate_,network,ipobody,xbody,ibody)
589 !
590  do jj=1,nope
591 !
592  j=jj
593 !
594  node1=kon(indexe+j)
595  jdof1=nactdof(0,node1)
596 !
597  do ll=jj,nope
598 !
599  l=ll
600 c m=0
601 !
602  node2=kon(indexe+l)
603  jdof2=nactdof(0,node2)
604 !
605 ! check whether one of the DOF belongs to a SPC or MPC
606 !
607  if((jdof1.gt.0).and.(jdof2.gt.0)) then
608  if(stiffonly(2).eq.1) then
609  call add_sm_st(au,ad,jq,irow,jdof1,jdof2,
610  & s(jj,ll),jj,ll)
611  else
612  call add_sm_ei(au,ad,aub,adb,jq,irow,jdof1,jdof2,
613  & s(jj,ll),sm(jj,ll),jj,ll)
614  endif
615  elseif((jdof1.gt.0).or.(jdof2.gt.0)) then
616 !
617 ! idof1: genuine DOF
618 ! idof2: nominal DOF of the SPC/MPC
619 !
620  if(jdof1.le.0) then
621  idof1=jdof2
622 c idof2=(node1-1)*8
623  idof2=jdof1
624  else
625  idof1=jdof1
626 c idof2=(node2-1)*8
627  idof2=jdof2
628  endif
629  if(nmpc.gt.0) then
630 c call nident(ikmpc,idof2,nmpc,id)
631 c if((id.gt.0).and.(ikmpc(id).eq.idof2)) then
632  if(idof2.ne.2*(idof2/2)) then
633 !
634 ! regular DOF / MPC
635 !
636 c id=ilmpc(id)
637  id=(-idof2+1)/2
638  ist=ipompc(id)
639  index=nodempc(3,ist)
640  if(index.eq.0) cycle
641  do
642  idof2=nactdof(nodempc(2,index),nodempc(1,index))
643  value=-coefmpc(index)*s(jj,ll)/coefmpc(ist)
644  if(idof1.eq.idof2) value=2.d0*value
645  if(idof2.gt.0) then
646  if(stiffonly(2).eq.1) then
647  call add_sm_st(au,ad,jq,irow,idof1,
648  & idof2,value,i0,i0)
649  else
650  valu2=-coefmpc(index)*sm(jj,ll)/
651  & coefmpc(ist)
652 c
653  if(idof1.eq.idof2) valu2=2.d0*valu2
654 c
655  call add_sm_ei(au,ad,aub,adb,jq,irow,
656  & idof1,idof2,value,valu2,i0,i0)
657  endif
658  endif
659  index=nodempc(3,index)
660  if(index.eq.0) exit
661  enddo
662  cycle
663  endif
664  endif
665 !
666 ! regular DOF / SPC
667 !
668  if(rhsi.eq.1) then
669  elseif(nmethod.eq.2) then
670  value=s(jj,ll)
671 c call nident(ikboun,idof2,nboun,id)
672 c icolumn=neq(2)+ilboun(id)
673  icolumn=neq(2)-idof2/2
674  call add_bo_st(au,jq,irow,idof1,icolumn,value)
675  endif
676  else
677 c idof1=(node1-1)*8
678 c idof2=(node2-1)*8
679  idof1=jdof1
680  idof2=jdof2
681  mpc1=0
682  mpc2=0
683  if(nmpc.gt.0) then
684 c call nident(ikmpc,idof1,nmpc,id1)
685 c if((id1.gt.0).and.(ikmpc(id1).eq.idof1)) mpc1=1
686 c call nident(ikmpc,idof2,nmpc,id2)
687 c if((id2.gt.0).and.(ikmpc(id2).eq.idof2)) mpc2=1
688  if(idof1.ne.2*(idof1/2)) mpc1=1
689  if(idof2.ne.2*(idof2/2)) mpc2=1
690  endif
691  if((mpc1.eq.1).and.(mpc2.eq.1)) then
692 c id1=ilmpc(id1)
693 c id2=ilmpc(id2)
694  id1=(-idof1+1)/2
695  id2=(-idof2+1)/2
696  if(id1.eq.id2) then
697 !
698 ! MPC id1 / MPC id1
699 !
700  ist=ipompc(id1)
701  index1=nodempc(3,ist)
702  if(index1.eq.0) cycle
703  do
704  idof1=nactdof(nodempc(2,index1),
705  & nodempc(1,index1))
706  index2=index1
707  do
708  idof2=nactdof(nodempc(2,index2),
709  & nodempc(1,index2))
710  value=coefmpc(index1)*coefmpc(index2)*
711  & s(jj,ll)/coefmpc(ist)/coefmpc(ist)
712  if((idof1.gt.0).and.(idof2.gt.0)) then
713  if(stiffonly(2).eq.1) then
714  call add_sm_st(au,ad,jq,irow,
715  & idof1,idof2,value,i0,i0)
716  else
717  valu2=coefmpc(index1)*coefmpc(index2)*
718  & sm(jj,ll)/coefmpc(ist)/coefmpc(ist)
719  call add_sm_ei(au,ad,aub,adb,jq,
720  & irow,idof1,idof2,value,valu2,i0,i0)
721  endif
722  endif
723 !
724  index2=nodempc(3,index2)
725  if(index2.eq.0) exit
726  enddo
727  index1=nodempc(3,index1)
728  if(index1.eq.0) exit
729  enddo
730  else
731 !
732 ! MPC id1 / MPC id2
733 !
734  ist1=ipompc(id1)
735  index1=nodempc(3,ist1)
736  if(index1.eq.0) cycle
737  do
738  idof1=nactdof(nodempc(2,index1),
739  & nodempc(1,index1))
740  ist2=ipompc(id2)
741  index2=nodempc(3,ist2)
742  if(index2.eq.0) then
743  index1=nodempc(3,index1)
744  if(index1.eq.0) then
745  exit
746  else
747  cycle
748  endif
749  endif
750  do
751  idof2=nactdof(nodempc(2,index2),
752  & nodempc(1,index2))
753  value=coefmpc(index1)*coefmpc(index2)*
754  & s(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
755  if(idof1.eq.idof2) value=2.d0*value
756  if((idof1.gt.0).and.(idof2.gt.0)) then
757  if(stiffonly(2).eq.1) then
758  call add_sm_st(au,ad,jq,irow,
759  & idof1,idof2,value,i0,i0)
760  else
761  valu2=coefmpc(index1)*coefmpc(index2)*
762  & sm(jj,ll)/coefmpc(ist1)/coefmpc(ist2)
763 c
764  if(idof1.eq.idof2) valu2=2.d0*valu2
765 c
766  call add_sm_ei(au,ad,aub,adb,jq,
767  & irow,idof1,idof2,value,valu2,i0,i0)
768  endif
769  endif
770 !
771  index2=nodempc(3,index2)
772  if(index2.eq.0) exit
773  enddo
774  index1=nodempc(3,index1)
775  if(index1.eq.0) exit
776  enddo
777  endif
778 c elseif(((mpc1.eq.1).or.(mpc2.eq.1)).and.rhsi)
779 c & then
780 c if(mpc1.eq.1) then
781 c!
782 c! MPC id1 / SPC
783 c!
784 c call nident(ikboun,idof2,nboun,id2)
785 c idof2=ilboun(id2)
786 c ist1=ipompc(id1)
787 c index1=nodempc(3,ist1)
788 c if(index1.eq.0) cycle
789 c do
790 c idof1=nactdof(nodempc(2,index1),
791 c & nodempc(1,index1))
792 c index1=nodempc(3,index1)
793 c if(index1.eq.0) exit
794 c enddo
795 c elseif(mpc2.eq.1) then
796 c!
797 c! MPC id2 / SPC
798 c!
799 c call nident(ikboun,idof1,nboun,id1)
800 c idof1=ilboun(id1)
801 c ist2=ipompc(id2)
802 c index2=nodempc(3,ist2)
803 c if(index2.eq.0) cycle
804 c endif
805  endif
806  endif
807  enddo
808 !
809  if(rhsi.eq.1) then
810 !
811 ! distributed forces
812 !
813  if(idist.ne.0) then
814  if(jdof1.le.0) then
815  if(nmpc.ne.0) then
816 c idof1=(node1-1)*8
817  idof1=jdof1
818 c call nident(ikmpc,idof1,nmpc,id)
819 c if((id.gt.0).and.(ikmpc(id).eq.idof1)) then
820  if(idof1.ne.2*(idof1/2)) then
821 c id=ilmpc(id)
822  id=(-idof1+1)/2
823  ist=ipompc(id)
824  index=nodempc(3,ist)
825  if(index.eq.0) cycle
826  do
827  jdof1=nactdof(nodempc(2,index),
828  & nodempc(1,index))
829  if(jdof1.gt.0) then
830  fext(jdof1)=fext(jdof1)
831  & -coefmpc(index)*ff(jj)
832  & /coefmpc(ist)
833  endif
834  index=nodempc(3,index)
835  if(index.eq.0) exit
836  enddo
837  endif
838  endif
839  cycle
840  endif
841  fext(jdof1)=fext(jdof1)+ff(jj)
842  endif
843  endif
844 !
845  enddo
846  enddo
847 !
848  endif
849 !
850  if(rhsi.eq.1) then
851 !
852 ! point forces
853 !
854  do i=1,nforc
855  if(ndirforc(i).gt.3) cycle
856 !
857 ! updating the external force vector for dynamic
858 ! calculations
859 !
860  if(nmethod.eq.4) fnext(ndirforc(i),nodeforc(1,i))=
861  & fnext(ndirforc(i),nodeforc(1,i))+xforc(i)
862 !
863  jdof=nactdof(ndirforc(i),nodeforc(1,i))
864  if(jdof.gt.0) then
865  fext(jdof)=fext(jdof)+xforc(i)
866  else
867 !
868 ! node is a dependent node of a MPC: distribute
869 ! the forces among the independent nodes
870 ! (proportional to their coefficients)
871 !
872 c jdof=8*(nodeforc(1,i)-1)+ndirforc(i)
873 c call nident(ikmpc,jdof,nmpc,id)
874 c if(id.gt.0) then
875 c if(ikmpc(id).eq.jdof) then
876  if(jdof.ne.2*(jdof/2)) then
877 c id=ilmpc(id)
878  id=(-jdof+1)/2
879  ist=ipompc(id)
880  index=nodempc(3,ist)
881  if(index.eq.0) cycle
882  do
883  jdof=nactdof(nodempc(2,index),nodempc(1,index))
884  if(jdof.gt.0) then
885  fext(jdof)=fext(jdof)-
886  & coefmpc(index)*xforc(i)/coefmpc(ist)
887  endif
888  index=nodempc(3,index)
889  if(index.eq.0) exit
890  enddo
891 c endif
892  endif
893  endif
894  enddo
895 !
896  endif
897 !
898  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
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)