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

Go to the source code of this file.

Functions/Subroutines

subroutine equationfs (inpc, textpart, ipompc, nodempc, coefmpc, nmpc, nmpc_, mpcfree, co, trab, ntrans, ikmpc, ilmpc, labmpc, istep, istat, n, iline, ipol, inl, ipoinp, inp, ipoinpc, lakon, ne, nload, sideload, ipkon, kon, nelemload)
 

Function/Subroutine Documentation

◆ equationfs()

subroutine equationfs ( character*1, dimension(*)  inpc,
character*132, dimension(16)  textpart,
integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  nmpc,
integer  nmpc_,
integer  mpcfree,
real*8, dimension(3,*)  co,
real*8, dimension(7,*)  trab,
integer  ntrans,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
character*20, dimension(*)  labmpc,
integer  istep,
integer  istat,
integer  n,
integer  iline,
integer  ipol,
integer  inl,
integer, dimension(2,*)  ipoinp,
integer, dimension(3,*)  inp,
integer, dimension(0:*)  ipoinpc,
character*8, dimension(*)  lakon,
integer  ne,
integer  nload,
character*20, dimension(*)  sideload,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
integer, dimension(2,*)  nelemload 
)
23 !
24 ! reading the input deck: *EQUATIONF
25 !
26  implicit none
27 !
28  character*1 inpc(*)
29  character*8 lakon(*)
30  character*20 labmpc(*),label,sideload(*)
31  character*132 textpart(16)
32 !
33  integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,istep,istat,
34  & n,i,j,ii,key,nterm,number,ntrans,ndir,indexe,
35  & mpcfreeold,ikmpc(*),ilmpc(*),id,idof,itr,iline,ipol,inl,
36  & ipoinp(2,*),inp(3,*),ipoinpc(0:*),
37  & k,m,iface,ifacel,nelem,ifaceq(8,6),ifacet(6,4),ifacew(8,5),
38  & loadid,ne,nope,nopes,ipkon(*),nload,nelemload(2,*),kon(*)
39 !
40  real*8 coefmpc(*),co(3,*),trab(7,*),a(3,3),x,cg(3)
41 !
42  data ifaceq /4,3,2,1,11,10,9,12,
43  & 5,6,7,8,13,14,15,16,
44  & 1,2,6,5,9,18,13,17,
45  & 2,3,7,6,10,19,14,18,
46  & 3,4,8,7,11,20,15,19,
47  & 4,1,5,8,12,17,16,20/
48  data ifacet /1,3,2,7,6,5,
49  & 1,2,4,5,9,8,
50  & 2,3,4,6,10,9,
51  & 1,4,3,8,10,7/
52  data ifacew /1,3,2,9,8,7,0,0,
53  & 4,5,6,10,11,12,0,0,
54  & 1,2,5,4,7,14,10,13,
55  & 2,3,6,5,8,15,11,14,
56  & 4,6,3,1,12,15,9,13/
57 !
58  do m=2,n
59  write(*,*)
60  & '*WARNING reading *EQUATIONF: parameter not recognized:'
61  write(*,*) ' ',
62  & textpart(m)(1:index(textpart(m),' ')-1)
63  call inputwarning(inpc,ipoinpc,iline,
64  &"*EQUATIONF%")
65  enddo
66 !
67  if(istep.gt.0) then
68  write(*,*)
69  & '*ERROR reading *EQUATIONF: *EQUATIONF should be placed'
70  write(*,*) ' before all step definitions'
71  call exit(201)
72  endif
73 !
74  do
75  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
76  & ipoinp,inp,ipoinpc)
77  if((istat.lt.0).or.(key.eq.1)) return
78  read(textpart(1)(1:10),'(i10)',iostat=istat) nterm
79 !
80  nmpc=nmpc+1
81  if(nmpc.gt.nmpc_) then
82  write(*,*) '*ERROR reading *EQUATIONF: increase nmpc_'
83  call exit(201)
84  endif
85 !
86  labmpc(nmpc)='FLUID '
87  ipompc(nmpc)=mpcfree
88  ii=0
89 !
90  do
91  call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl,
92  & ipoinp,inp,ipoinpc)
93  if((istat.lt.0).or.(key.eq.1)) then
94  write(*,*) '*ERROR reading *EQUATIONF: mpc definition ',
95  & nmpc
96  write(*,*) ' is not complete. '
97  call inputerror(inpc,ipoinpc,iline,
98  &"*EQUATIONF%")
99  call exit(201)
100  endif
101 !
102  do i=1,n/4
103 !
104  read(textpart((i-1)*4+1)(1:10),'(i10)',iostat=istat)nelem
105  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
106  &"*EQUATIONF%")
107  if((nelem.gt.ne).or.(nelem.le.0)) then
108  write(*,*) '*ERROR reading *EQUATIONF:'
109  write(*,*) ' element ',nelem,' is not defined'
110  call exit(201)
111  endif
112 !
113  read(textpart((i-1)*4+2)(2:2),'(i1)',iostat=istat)
114  & ifacel
115  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
116  &"*EQUATIONF%")
117  iface=10*nelem+ifacel
118 !
119  read(textpart((i-1)*4+3)(1:10),'(i10)',iostat=istat) ndir
120  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
121  &"*EQUATIONF%")
122  if(ndir.le.3) then
123  elseif(ndir.eq.4) then
124  write(*,*) '*ERROR in equationfs: an equation'
125  write(*,*) ' on DOF 4 is not allowed'
126  call exit(201)
127  elseif(ndir.eq.5) then
128  write(*,*) '*ERROR in equationfs: an equation'
129  write(*,*) ' on DOF 5 is not allowed'
130  call exit(201)
131  elseif(ndir.eq.6) then
132  write(*,*) '*ERROR in equationfs: an equation'
133  write(*,*) ' on DOF 6 is not allowed'
134  call exit(201)
135  elseif(ndir.eq.8) then
136  ndir=4
137  elseif(ndir.eq.11) then
138  ndir=0
139  else
140  write(*,*) '*ERROR reading *EQUATIONF:'
141  write(*,*) ' direction',ndir,' is not defined'
142  call exit(201)
143  endif
144 !
145  read(textpart((i-1)*4+4)(1:20),'(f20.0)',iostat=istat) x
146  if(istat.gt.0) call inputerror(inpc,ipoinpc,iline,
147  &"*EQUATIONF%")
148 !
149 ! check whether the face is transformed
150 !
151  if(ntrans.le.0) then
152  itr=0
153  else
154  label(1:20)='T '
155  write(label(2:2),'(i1)') ifacel
156  call identifytransform(nelem,label,nelemload,sideload,
157  & nload,loadid)
158  if(loadid.eq.0) then
159  itr=0
160  else
161  itr=nelemload(2,loadid)
162  endif
163  endif
164 !
165  if((itr.eq.0).or.(ndir.eq.0).or.(ndir.eq.4)) then
166  nodempc(1,mpcfree)=iface
167  nodempc(2,mpcfree)=ndir
168  coefmpc(mpcfree)=x
169 !
170 ! updating ikmpc and ilmpc
171 !
172  if(ii.eq.0) then
173  idof=-(8*(iface-1)+ndir)
174  call nident(ikmpc,idof,nmpc-1,id)
175  if(id.gt.0) then
176  if(ikmpc(id).eq.idof) then
177  write(*,100)
178  & (ikmpc(id))/8+1,ikmpc(id)-8*((ikmpc(id))/8)
179  call exit(201)
180  endif
181  endif
182  do j=nmpc,id+2,-1
183  ikmpc(j)=ikmpc(j-1)
184  ilmpc(j)=ilmpc(j-1)
185  enddo
186  ikmpc(id+1)=idof
187  ilmpc(id+1)=nmpc
188  endif
189 !
190  mpcfreeold=mpcfree
191  mpcfree=nodempc(3,mpcfree)
192  if(mpcfree.eq.0) then
193  write(*,*)
194  & '*ERROR reading *EQUATIONF: increase memmpc_'
195  call exit(201)
196  endif
197  else
198 !
199 ! transformation applies
200 !
201 ! number of nodes belonging to the face
202 !
203  if(lakon(nelem)(4:4).eq.'8') then
204  nope=8
205  nopes=4
206  elseif(lakon(nelem)(4:4).eq.'4') then
207  nope=4
208  nopes=3
209  elseif(lakon(nelem)(4:4).eq.'6') then
210  nope=6
211  if(ifacel.le.2) then
212  nopes=3
213  else
214  nopes=4
215  endif
216  endif
217 !
218 ! determining the center of gravity
219 !
220  do j=1,3
221  cg(j)=0.d0
222  enddo
223 !
224  indexe=ipkon(nelem)
225  if(nope.eq.8) then
226  do k=1,nopes
227  do j=1,3
228  cg(j)=cg(j)+
229  & co(j,kon(indexe+ifaceq(k,ifacel)))
230  enddo
231  enddo
232  elseif(nope.eq.4) then
233  do k=1,nopes
234  do j=1,3
235  cg(j)=cg(j)+
236  & co(j,kon(indexe+ifacet(k,ifacel)))
237  enddo
238  enddo
239  else
240  do k=1,nopes
241  do j=1,3
242  cg(j)=cg(j)+
243  & co(j,kon(indexe+ifacew(k,ifacel)))
244  enddo
245  enddo
246  endif
247  do j=1,3
248  cg(j)=cg(j)/nopes
249  enddo
250 !
251 ! determining the transformation coefficients at the center
252 ! of gravity
253 !
254  call transformatrix(trab(1,itr),
255  & cg,a)
256 !
257  number=ndir-1
258  if(ii.eq.0) then
259 !
260 ! determining which direction to use for the
261 ! dependent side: should not occur on the dependent
262 ! side in another MPC and should have a nonzero
263 ! coefficient
264 !
265  do j=1,3
266  number=number+1
267  if(number.gt.3) number=1
268  idof=-(8*(iface-1)+number)
269  call nident(ikmpc,idof,nmpc-1,id)
270  if(id.gt.0) then
271  if(ikmpc(id).eq.idof) then
272  cycle
273  endif
274  endif
275  if(dabs(a(number,ndir)).lt.1.d-5) cycle
276  exit
277  enddo
278  if(j.gt.3) then
279  write(*,*)
280  & '*ERROR reading *EQUATIONF: MPC on face'
281  write(*,*) ifacel,' of element',nelem
282  write(*,*) ' in transformed coordinates'
283  write(*,*) ' cannot be converted in MPC: all'
284  write(*,*) ' DOFs in the node are used as'
285  write(*,*) ' dependent nodes in other MPCs'
286  call exit(201)
287  endif
288  number=number-1
289 !
290 ! updating ikmpc and ilmpc
291 !
292  do j=nmpc,id+2,-1
293  ikmpc(j)=ikmpc(j-1)
294  ilmpc(j)=ilmpc(j-1)
295  enddo
296  ikmpc(id+1)=idof
297  ilmpc(id+1)=nmpc
298  endif
299 !
300  do j=1,3
301  number=number+1
302  if(number.gt.3) number=1
303  if(dabs(a(number,ndir)).lt.1.d-5) cycle
304  nodempc(1,mpcfree)=iface
305  nodempc(2,mpcfree)=number
306  coefmpc(mpcfree)=x*a(number,ndir)
307  mpcfreeold=mpcfree
308  mpcfree=nodempc(3,mpcfree)
309  if(mpcfree.eq.0) then
310  write(*,*)
311  & '*ERROR reading *EQUATIONF: increase memmpc_'
312  call exit(201)
313  endif
314  enddo
315  endif
316 !
317  ii=ii+1
318  enddo
319 !
320  if(ii.eq.nterm) then
321  nodempc(3,mpcfreeold)=0
322  exit
323  endif
324  enddo
325  enddo
326 !
327  100 format(/,'*ERROR reading *EQUATIONF: the DOF corresponding to',
328  & /,'iface ',i1,' of element',i10,' in direction',
329  & i5,' is detected on',
330  & /,'the dependent side of two different MPC''s')
331  return
inputwarning
subroutine inputwarning(inpc, ipoinpc, iline, text)
Definition: inputwarning.f:20
identifytransform
subroutine identifytransform(nelement, label, nelemload, sideload, nload, loadid)
Definition: identifytransform.f:21
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
getnewline
subroutine getnewline(inpc, textpart, istat, n, key, iline, ipol, inl, ipoinp, inp, ipoinpc)
Definition: getnewline.f:21
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
inputerror
subroutine inputerror(inpc, ipoinpc, iline, text)
Definition: inputerror.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)