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

Go to the source code of this file.

Functions/Subroutines

subroutine knotmpc (ipompc, nodempc, coefmpc, irefnode, irotnode, iexpnode, labmpc, nmpc, nmpc_, mpcfree, ikmpc, ilmpc, nk, nk_, nodeboun, ndirboun, ikboun, ilboun, nboun, nboun_, idepnodes, typeboun, co, xboun, istep, k, ndepnodes, idim, e1, e2, t1)
 

Function/Subroutine Documentation

◆ knotmpc()

subroutine knotmpc ( integer, dimension(*)  ipompc,
integer, dimension(3,*)  nodempc,
real*8, dimension(*)  coefmpc,
integer  irefnode,
integer  irotnode,
integer  iexpnode,
character*20, dimension(*)  labmpc,
integer  nmpc,
integer  nmpc_,
integer  mpcfree,
integer, dimension(*)  ikmpc,
integer, dimension(*)  ilmpc,
integer  nk,
integer  nk_,
integer, dimension(*)  nodeboun,
integer, dimension(*)  ndirboun,
integer, dimension(*)  ikboun,
integer, dimension(*)  ilboun,
integer  nboun,
integer  nboun_,
integer, dimension(*)  idepnodes,
character*1, dimension(*)  typeboun,
real*8, dimension(3,*)  co,
real*8, dimension(*)  xboun,
integer  istep,
integer  k,
integer  ndepnodes,
integer  idim,
real*8, dimension(3)  e1,
real*8, dimension(3)  e2,
real*8, dimension(3)  t1 
)
24 !
25 ! generates three knot MPC's for node "node" about reference
26 ! (translational) node irefnode and rotational node irotnode
27 !
28  implicit none
29 !
30  character*1 typeboun(*)
31  character*20 labmpc(*)
32 !
33  integer ipompc(*),nodempc(3,*),nmpc,nmpc_,mpcfree,nk,nk_,
34  & ikmpc(*),idepnodes(*),k,ndepnodes,idim,n,matz,ier,i,
35  & ilmpc(*),node,id,mpcfreeold,j,idof,l,nodeboun(*),
36  & ndirboun(*),ikboun(*),ilboun(*),nboun,nboun_,irefnode,
37  & irotnode,iexpnode,istep,ispcnode,inode
38 !
39  real*8 coefmpc(*),co(3,*),xboun(*),e(3,3,3),dc(3,3,3),s(3,3),
40  & w(3),z(3,3),fv1(3),fv2(3),e1(3),e2(3),t1(3),sx,sy,sz,sxx,
41  & sxy,sxz,syy,syz,szz,u2(3,3),u3(3,3)
42 !
43 ! e_ijk symbol
44 !
45  data e /0.,0.,0.,0.,0.,-1.,0.,1.,0.,
46  & 0.,0.,1.,0.,0.,0.,-1.,0.,0.,
47  & 0.,-1.,0.,1.,0.,0.,0.,0.,0./
48 !
49 ! dc_ijk=e_ikj
50 !
51  data dc /0.,0.,0.,0.,0.,1.,0.,-1.,0.,
52  & 0.,0.,-1.,0.,0.,0.,1.,0.,0.,
53  & 0.,1.,0.,-1.,0.,0.,0.,0.,0./
54 !
55 c idim=1
56 !
57 ! on entry: if idim=1: only 2-d elements
58 ! if idim=3: at least one 1-d element
59 !
60  if((istep.gt.1).and.(k.eq.1)) then
61 !
62 ! determine the dimensionality of the knot
63 ! (only for step 2 or higher since in step 1 the geometry
64 ! is not known yet)
65 !
66 ! determine the area moments of inertia
67 !
68  sx=0.d0
69  sy=0.d0
70  sz=0.d0
71  sxx=0.d0
72  sxy=0.d0
73  sxz=0.d0
74  syy=0.d0
75  syz=0.d0
76  szz=0.d0
77 !
78  do i=1,ndepnodes
79  node=idepnodes(i)
80  sx=sx+co(1,node)
81  sy=sy+co(2,node)
82  sz=sz+co(3,node)
83  sxx=sxx+co(1,node)*co(1,node)
84  sxy=sxy+co(1,node)*co(2,node)
85  sxz=sxz+co(1,node)*co(3,node)
86  syy=syy+co(2,node)*co(2,node)
87  syz=syz+co(2,node)*co(3,node)
88  szz=szz+co(3,node)*co(3,node)
89  enddo
90 !
91  sxx=sxx-sx*sx/ndepnodes
92  sxy=sxy-sx*sy/ndepnodes
93  sxz=sxz-sx*sz/ndepnodes
94  syy=syy-sy*sy/ndepnodes
95  syz=syz-sy*sz/ndepnodes
96  szz=szz-sz*sz/ndepnodes
97 !
98 c write(*,*) 'sxx...',sxx,sxy,sxz,syy,syz,szz
99 !
100  s(1,1)=sxx
101  s(1,2)=sxy
102  s(1,3)=sxz
103  s(2,1)=sxy
104  s(2,2)=syy
105  s(2,3)=syz
106  s(3,1)=sxz
107  s(3,2)=syz
108  s(3,3)=szz
109 !
110 ! determining the eigenvalues
111 !
112  n=3
113  matz=1
114  ier=0
115  call rs(n,n,s,w,matz,z,fv1,fv2,ier)
116  if(ier.ne.0) then
117  write(*,*) '*ERROR in knotmpc while calculating the'
118  write(*,*) ' eigenvalues/eigenvectors'
119  call exit(201)
120  endif
121 !
122 ! the eigenvalues are the moments of inertia w.r.t. the
123 ! plane orthogonal to the eigenvector
124 !
125 ! dimension=1 if the two lowest eigenvalues are zero
126 ! dimension=2 if only the lowest eigenvalue is zero
127 ! else dimension=3
128 !
129 c write(*,*) 'eigenvalues ',w(1),w(2),w(3)
130  if((w(1).lt.1.d-10).and.(w(2).lt.1.d-10)) then
131  idim=min(idim,1)
132 c idim=1
133  elseif(w(1).lt.1.d-10) then
134  idim=min(idim,2)
135 c idim=2
136  else
137  idim=min(idim,1)
138 c idim=3
139  endif
140 c write(*,*) 'knotmpc irefnode= ',irefnode,' idim= ',idim
141 !
142 ! defining a local coordinate system for idim=2
143 !
144  if(idim.eq.2) then
145  do i=1,3
146  t1(i)=z(i,1)
147  e2(i)=z(i,2)
148  e1(i)=z(i,3)
149  enddo
150 !
151 ! check whether e1-e2-t1 is a rhs system
152 !
153  if(t1(1)*(e1(2)*e2(3)-e1(3)*e2(2))-
154  & t1(2)*(e1(1)*e2(3)-e1(3)*e2(1))+
155  & t1(3)*(e1(1)*e2(2)-e1(2)*e2(1)).lt.0.d0) then
156  do i=1,3
157  t1(i)=-t1(i)
158  enddo
159  endif
160 c write(*,*) 't1 ',t1(1),t1(2),t1(3)
161 c write(*,*) 'e1 ',e1(1),e1(2),e1(3)
162 c write(*,*) 'e2 ',e2(1),e2(2),e2(3)
163 !
164 ! storing t1 and e1 as coordinates of irotnode and
165 ! iexpnode, respectively
166 !
167  do i=1,3
168  co(i,irotnode)=t1(i)
169  co(i,iexpnode)=e1(i)
170  enddo
171  endif
172  endif
173 !
174  inode=idepnodes(k)
175 !
176  nk=nk+1
177  if(nk.gt.nk_) then
178  write(*,*) '*ERROR in knotmpc: increase nk_'
179  call exit(201)
180  endif
181 !
182  ispcnode=nk
183 !
184  if((idim.eq.1).or.(idim.eq.3)) then
185 !
186 ! knot nodes lie on a line
187 !
188  do j=1,3
189  idof=8*(inode-1)+j
190  call nident(ikmpc,idof,nmpc,id)
191  if(id.gt.0) then
192  if(ikmpc(id).eq.idof) then
193  cycle
194  endif
195  endif
196  nmpc=nmpc+1
197  if(nmpc.gt.nmpc_) then
198  write(*,*) '*ERROR in knotmpc: increase nmpc_'
199  call exit(201)
200  endif
201 !
202  ipompc(nmpc)=mpcfree
203  labmpc(nmpc)='KNOT '
204  write(labmpc(nmpc)(5:5),'(i1)') idim
205 !
206  do l=nmpc,id+2,-1
207  ikmpc(l)=ikmpc(l-1)
208  ilmpc(l)=ilmpc(l-1)
209  enddo
210  ikmpc(id+1)=idof
211  ilmpc(id+1)=nmpc
212 !
213  nodempc(1,mpcfree)=inode
214  nodempc(2,mpcfree)=j
215  coefmpc(mpcfree)=1.d0
216  mpcfree=nodempc(3,mpcfree)
217 !
218 ! translation term
219 !
220  nodempc(1,mpcfree)=irefnode
221  nodempc(2,mpcfree)=j
222  coefmpc(mpcfree)=-1.d0
223  mpcfree=nodempc(3,mpcfree)
224 !
225 ! expansion term
226 !
227  nodempc(1,mpcfree)=iexpnode
228  nodempc(2,mpcfree)=1
229  if(istep.gt.1) then
230  coefmpc(mpcfree)=co(j,irefnode)-co(j,inode)
231  endif
232  mpcfree=nodempc(3,mpcfree)
233 !
234 ! rotation terms
235 !
236  nodempc(1,mpcfree)=irotnode
237  nodempc(2,mpcfree)=1
238  if(istep.gt.1) then
239  coefmpc(mpcfree)=dc(j,1,1)*(co(1,irefnode)-co(1,inode))+
240  & dc(j,2,1)*(co(2,irefnode)-co(2,inode))+
241  & dc(j,3,1)*(co(3,irefnode)-co(3,inode))
242  endif
243  mpcfree=nodempc(3,mpcfree)
244  nodempc(1,mpcfree)=irotnode
245  nodempc(2,mpcfree)=2
246  if(istep.gt.1) then
247  coefmpc(mpcfree)=dc(j,1,2)*(co(1,irefnode)-co(1,inode))+
248  & dc(j,2,2)*(co(2,irefnode)-co(2,inode))+
249  & dc(j,3,2)*(co(3,irefnode)-co(3,inode))
250  endif
251  mpcfree=nodempc(3,mpcfree)
252  nodempc(1,mpcfree)=irotnode
253  nodempc(2,mpcfree)=3
254  if(istep.gt.1) then
255  coefmpc(mpcfree)=dc(j,1,3)*(co(1,irefnode)-co(1,inode))+
256  & dc(j,2,3)*(co(2,irefnode)-co(2,inode))+
257  & dc(j,3,3)*(co(3,irefnode)-co(3,inode))
258  endif
259  mpcfree=nodempc(3,mpcfree)
260  nodempc(1,mpcfree)=ispcnode
261  nodempc(2,mpcfree)=j
262  coefmpc(mpcfree)=1.d0
263  mpcfreeold=mpcfree
264  mpcfree=nodempc(3,mpcfree)
265  nodempc(3,mpcfreeold)=0
266  idof=8*(ispcnode-1)+j
267  call nident(ikboun,idof,nboun,id)
268  nboun=nboun+1
269  if(nboun.gt.nboun_) then
270  write(*,*) '*ERROR in knotmpc: increase nboun_'
271  call exit(201)
272  endif
273  nodeboun(nboun)=ispcnode
274  ndirboun(nboun)=j
275  typeboun(nboun)='R'
276  if(istep.gt.1) then
277  xboun(nboun)=0.d0
278  endif
279  do l=nboun,id+2,-1
280  ikboun(l)=ikboun(l-1)
281  ilboun(l)=ilboun(l-1)
282  enddo
283  ikboun(id+1)=idof
284  ilboun(id+1)=nboun
285  enddo
286  elseif(idim.eq.2) then
287 !
288 ! knot nodes lie in a plane
289 !
290  do i=1,3
291  do j=1,3
292  u2(i,j)=2.d0*e1(i)*e1(j)
293  u3(i,j)=2.d0*e2(i)*e2(j)
294  enddo
295  enddo
296 !
297  do j=1,3
298  idof=8*(inode-1)+j
299  call nident(ikmpc,idof,nmpc,id)
300  if(id.gt.0) then
301  if(ikmpc(id).eq.idof) then
302  cycle
303  endif
304  endif
305  nmpc=nmpc+1
306  if(nmpc.gt.nmpc_) then
307  write(*,*) '*ERROR in knotmpc: increase nmpc_'
308  call exit(201)
309  endif
310 !
311  ipompc(nmpc)=mpcfree
312  labmpc(nmpc)='KNOT2 '
313 !
314  do l=nmpc,id+2,-1
315  ikmpc(l)=ikmpc(l-1)
316  ilmpc(l)=ilmpc(l-1)
317  enddo
318  ikmpc(id+1)=idof
319  ilmpc(id+1)=nmpc
320 !
321  nodempc(1,mpcfree)=inode
322  nodempc(2,mpcfree)=j
323  coefmpc(mpcfree)=1.d0
324  mpcfree=nodempc(3,mpcfree)
325 !
326 ! translation term
327 !
328  nodempc(1,mpcfree)=irefnode
329  nodempc(2,mpcfree)=j
330  coefmpc(mpcfree)=-1.d0
331  mpcfree=nodempc(3,mpcfree)
332 !
333 ! expansion terms
334 !
335 ! first term is "removed" (= amalgated with the second term)
336 ! since u1 is the null matrix
337 !
338  nodempc(1,mpcfree)=iexpnode
339  nodempc(2,mpcfree)=2
340  coefmpc(mpcfree)=0.d0
341  mpcfree=nodempc(3,mpcfree)
342 !
343  nodempc(1,mpcfree)=iexpnode
344  nodempc(2,mpcfree)=2
345  coefmpc(mpcfree)=u2(j,1)*(co(1,irefnode)-co(1,inode))+
346  & u2(j,2)*(co(2,irefnode)-co(2,inode))+
347  & u2(j,3)*(co(3,irefnode)-co(3,inode))
348  mpcfree=nodempc(3,mpcfree)
349 !
350  nodempc(1,mpcfree)=iexpnode
351  nodempc(2,mpcfree)=3
352  coefmpc(mpcfree)=u3(j,1)*(co(1,irefnode)-co(1,inode))+
353  & u3(j,2)*(co(2,irefnode)-co(2,inode))+
354  & u3(j,3)*(co(3,irefnode)-co(3,inode))
355  mpcfree=nodempc(3,mpcfree)
356 !
357 ! rotation terms
358 !
359  nodempc(1,mpcfree)=irotnode
360  nodempc(2,mpcfree)=1
361  if(istep.gt.1) then
362  coefmpc(mpcfree)=dc(j,1,1)*(co(1,irefnode)-co(1,inode))+
363  & dc(j,2,1)*(co(2,irefnode)-co(2,inode))+
364  & dc(j,3,1)*(co(3,irefnode)-co(3,inode))
365  endif
366  mpcfree=nodempc(3,mpcfree)
367  nodempc(1,mpcfree)=irotnode
368  nodempc(2,mpcfree)=2
369  if(istep.gt.1) then
370  coefmpc(mpcfree)=dc(j,1,2)*(co(1,irefnode)-co(1,inode))+
371  & dc(j,2,2)*(co(2,irefnode)-co(2,inode))+
372  & dc(j,3,2)*(co(3,irefnode)-co(3,inode))
373  endif
374  mpcfree=nodempc(3,mpcfree)
375  nodempc(1,mpcfree)=irotnode
376  nodempc(2,mpcfree)=3
377  if(istep.gt.1) then
378  coefmpc(mpcfree)=dc(j,1,3)*(co(1,irefnode)-co(1,inode))+
379  & dc(j,2,3)*(co(2,irefnode)-co(2,inode))+
380  & dc(j,3,3)*(co(3,irefnode)-co(3,inode))
381  endif
382  mpcfree=nodempc(3,mpcfree)
383  nodempc(1,mpcfree)=ispcnode
384  nodempc(2,mpcfree)=j
385  coefmpc(mpcfree)=1.d0
386  mpcfreeold=mpcfree
387  mpcfree=nodempc(3,mpcfree)
388  nodempc(3,mpcfreeold)=0
389  idof=8*(ispcnode-1)+j
390  call nident(ikboun,idof,nboun,id)
391  nboun=nboun+1
392  if(nboun.gt.nboun_) then
393  write(*,*) '*ERROR in knotmpc: increase nboun_'
394  call exit(201)
395  endif
396  nodeboun(nboun)=ispcnode
397  ndirboun(nboun)=j
398  typeboun(nboun)='R'
399  if(istep.gt.1) then
400  xboun(nboun)=0.d0
401  endif
402  do l=nboun,id+2,-1
403  ikboun(l)=ikboun(l-1)
404  ilboun(l)=ilboun(l-1)
405  enddo
406  ikboun(id+1)=idof
407  ilboun(id+1)=nboun
408  enddo
409  else
410 !
411 ! to do
412 !
413  endif
414 !
415  return
min
#define min(a, b)
Definition: cascade.c:31
rs
subroutine rs(nm, n, a, w, matz, z, fv1, fv2, ierr)
Definition: rs.f:27
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
Hosted by OpenAircraft.com, (Michigan UAV, LLC)