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

Go to the source code of this file.

Functions/Subroutines

subroutine normalsforequ_se (nk, co, iponoelfa, inoelfa, konfa, ipkonfa, lakonfa, ne, iponor, xnor, nodedesiinv, jobnamef)
 

Function/Subroutine Documentation

◆ normalsforequ_se()

subroutine normalsforequ_se ( integer  nk,
real*8, dimension(3,*)  co,
integer, dimension(*)  iponoelfa,
integer, dimension(3,*)  inoelfa,
integer, dimension(*)  konfa,
integer, dimension(*)  ipkonfa,
character*8, dimension(*)  lakonfa,
integer  ne,
integer, dimension(*)  iponor,
real*8, dimension(*)  xnor,
integer, dimension(*)  nodedesiinv,
character*132  jobnamef 
)
21 !
22 ! calculates normals on surface
23 !
24  implicit none
25 !
26  character*132 jobnamef,fnequ
27  character*8 lakonfa(*),lakonfaloc
28 !
29  integer nk,iponoelfa(*),inoelfa(3,*),konfa(*),ipkonfa(*),ne,
30  & i,ndepnodes,index,nexp,nnor,nel,ielem,indexe,j,iel(100),
31  & jl(100),ial(100),ifi(100),indexx,k,l,ifix,nemin,jact,ixfree,
32  & node,idim,iponor(*),nodedesiinv(*),len,ndet(3),nsort(3)
33 !
34  real*8 co(3,*),xnor(*),xno(3,100),xi,et,coloc6(2,6),coloc8(2,8),
35  & xl(3,8),dd,xnoref(3),dot,xnorloc(3,3),det(3),sort(3)
36 !
37  data coloc6 /0.d0,0.d0,1.d0,0.d0,0.d0,1.d0,0.5d0,0.d0,
38  & 0.5d0,0.5d0,0.d0,0.5d0/
39  data coloc8 /-1.d0,-1.d0,1.d0,-1.d0,1.d0,1.d0,-1.d0,1.d0,
40  & 0.d0,-1.d0,1.d0,0.d0,0.d0,1.d0,-1.d0,0.d0/
41 !
42 ! calculating the normals in nodes belonging to shells
43 !
44  do len=1,132
45  if(jobnamef(len:len).eq.' ') exit
46  enddo
47  len=len-1
48 
49  fnequ=jobnamef(1:len)//'.equ'
50  open(20,file=fnequ(1:len+4),status='unknown',err=100)
51  close(20,status='delete',err=101)
52  open(20,file=fnequ(1:len+4),status='unknown',err=100)
53  write(20,102)
54  write(20,103)
55  102 format('**SUMMARY OF EUQATIONS FOR MESH-UPDATE ')
56  103 format('*EQUATION')
57 !
58  ixfree=0
59 !
60  do i=1,nk
61  if(i.eq.1699) then
62  ndepnodes=0
63  endif
64  ndepnodes=0
65  idim=0
66  index=iponoelfa(i)
67  if(index.eq.0) cycle
68 !
69 ! nexp indicates how many times the node was expanded
70 !
71  nexp=0
72 !
73 ! nnor indicates whether the expanded nodes lie on a point
74 ! (nnor=0, only for plane stress, plane strain or axisymmetric
75 ! elements), on a line (nnor=1) or in a plane (nnor=2)
76 !
77  nnor=0
78 !
79 ! locating the shell elements to which node i belongs
80 !
81  nel=0
82  do
83  if(index.eq.0) exit
84  ielem=inoelfa(1,index)
85  if(lakonfa(ielem)(1:1).eq.'S') nnor=1
86  indexe=ipkonfa(ielem)
87  nel=nel+1
88  if(nel.gt.100) then
89  write(*,*) '*ERROR in normalsforequ_se: more '
90  write(*,*) ' than 100 shell elements '
91  write(*,*) ' share the same node'
92  call exit(201)
93  endif
94  j=inoelfa(2,index)
95  jl(nel)=j
96  iel(nel)=ielem
97  index=inoelfa(3,index)
98  enddo
99 !
100  if(nel.gt.0) then
101  do j=1,nel
102  ial(j)=0
103  enddo
104 !
105 ! estimate the normal
106 !
107  do j=1,nel
108  indexe=ipkonfa(iel(j))
109  indexx=iponor(indexe+jl(j))
110  if(indexx.ge.0) then
111  do k=1,3
112  xno(k,j)=xnor(indexx+k)
113  enddo
114  ifi(j)=1
115  cycle
116  else
117  ifi(j)=0
118  endif
119 !
120 ! local normal on the element (Jacobian)
121 !
122  lakonfaloc=lakonfa(iel(j))
123  if(lakonfa(iel(j))(2:2).eq.'3') then
124  xi=coloc6(1,jl(j))
125  et=coloc6(2,jl(j))
126  do k=1,3
127  node=konfa(indexe+k)
128  do l=1,3
129  xl(l,k)=co(l,node)
130  enddo
131  enddo
132  call norshell3(xi,et,xl,xno(1,j))
133  elseif(lakonfa(iel(j))(2:2).eq.'4') then
134  xi=coloc8(1,jl(j))
135  et=coloc8(2,jl(j))
136  do k=1,4
137  node=konfa(indexe+k)
138  do l=1,3
139  xl(l,k)=co(l,node)
140  enddo
141  enddo
142  call norshell4(xi,et,xl,xno(1,j))
143  elseif(lakonfa(iel(j))(2:2).eq.'6') then
144  xi=coloc6(1,jl(j))
145  et=coloc6(2,jl(j))
146  do k=1,6
147  node=konfa(indexe+k)
148  do l=1,3
149  xl(l,k)=co(l,node)
150  enddo
151  enddo
152  call norshell6(xi,et,xl,xno(1,j))
153  elseif(lakonfa(iel(j))(2:2).eq.'8') then
154  xi=coloc8(1,jl(j))
155  et=coloc8(2,jl(j))
156  do k=1,8
157  node=konfa(indexe+k)
158  do l=1,3
159  xl(l,k)=co(l,node)
160  enddo
161  enddo
162  call norshell8(xi,et,xl,xno(1,j))
163  endif
164 !
165  dd=dsqrt(xno(1,j)**2+xno(2,j)**2+xno(3,j)**2)
166  if(dd.lt.1.d-10) then
167  write(*,*) '*ERROR in normalsforequ_se: size '
168  write(*,*) ' of estimatedshell normal in
169  &node ',i,' element ',iel(j)
170  write(*,*) ' is smaller than 1.e-10'
171  call exit(201)
172  endif
173  do k=1,3
174  xno(k,j)=xno(k,j)/dd
175  enddo
176  enddo
177 !
178  do
179 !
180 ! determining a fixed normal which was not treated yet,
181 ! or, if none is left, the minimum element number of all
182 ! elements containing node i and for which no normal was
183 ! determined yet
184 !
185 ! if ial(j)=0: the normal on this element has not been
186 ! treated yet
187 ! if ial(j)=2: normal has been treated
188 !
189  ifix=0
190  nemin=ne+1
191  do j=1,nel
192  if(ial(j).ne.0) cycle
193  if(ifi(j).eq.1) then
194  jact=j
195  ifix=1
196  exit
197  endif
198  enddo
199  if(ifix.eq.0) then
200  do j=1,nel
201  if(ial(j).eq.0) then
202  if(iel(j).lt.nemin) then
203  nemin=iel(j)
204  jact=j
205  endif
206  endif
207  enddo
208  if(nemin.eq.ne+1) exit
209  endif
210 !
211  do j=1,3
212  xnoref(j)=xno(j,jact)
213  enddo
214 !
215 ! determining all elements whose normal in node i makes an
216 ! angle smaller than 0.5 or 20 degrees with the reference normal,
217 ! depending whether the reference normal was given by the
218 ! user or is being calculated; the thickness and offset must
219 ! also fit.
220 !
221 ! if ial(j)=1: normal on element is being treated now
222 !
223  do j=1,nel
224  if(ial(j).eq.2) cycle
225  if(j.eq.jact) then
226  ial(jact)=1
227  else
228  dot=xno(1,j)*xnoref(1)+xno(2,j)*xnoref(2)+
229  & xno(3,j)*xnoref(3)
230  if(ifix.eq.0) then
231  if(dot.gt.0.939693d0)then
232  if((lakonfa(iel(j))(1:3).eq.
233  & lakonfa(iel(jact))(1:3))
234  & .or.
235  & ((lakonfa(iel(j))(1:1).eq.'S').and.
236  & (lakonfa(iel(jact))(1:1).eq.'S')))
237  & ial(j)=1
238 c
239  if(dot.lt.0.999962d0) nnor=2
240 c
241  else
242  if((lakonfa(iel(j))(1:1).eq.'S').and.
243  & (lakonfa(iel(jact))(1:1).eq.'S')) then
244 !
245 ! if the normals have the opposite
246 ! direction, the expanded nodes are on a
247 ! straight line
248 !
249  if(dot.gt.-0.999962) then
250  nnor=2
251  else
252  write(*,*) '*INFO in gen3dnor: in some
253  &nodes opposite normals are defined'
254  endif
255  endif
256  endif
257  else
258  if(dot.gt.0.999962d0) then
259  if((lakonfa(iel(j))(1:3).eq.
260  & lakonfa(iel(jact))(1:3))
261  & .or.
262  & ((lakonfa(iel(j))(1:1).eq.'S').and.
263  & (lakonfa(iel(jact))(1:1).eq.'S')))
264  & ial(j)=1
265  else
266  if((lakonfa(iel(j))(1:1).eq.'S').and.
267  & (lakonfa(iel(jact))(1:1).eq.'S')) then
268 !
269 ! if the normals have the opposite
270 ! direction, the expanded nodes are on a
271 ! straight line
272 !
273  if(dot.gt.-0.999962) then
274  nnor=2
275  else
276  write(*,*) '*INFO in gen3dnor: in some
277  &nodes opposite normals are defined'
278  endif
279  endif
280  endif
281  endif
282  endif
283  enddo
284 !
285 ! determining the mean normal for the selected elements
286 !
287  if(ifix.eq.0) then
288  do j=1,3
289  xnoref(j)=0.d0
290  enddo
291  do j=1,nel
292  if(ial(j).eq.1) then
293  do k=1,3
294  xnoref(k)=xnoref(k)+xno(k,j)
295  enddo
296  endif
297  enddo
298  dd=dsqrt(xnoref(1)**2+xnoref(2)**2+xnoref(3)**2)
299  if(dd.lt.1.d-10) then
300  write(*,*) '*ERROR in gen3dnor: size of'
301  write(*,*) ' estimated shell normal is'
302  write(*,*) ' smaller than 1.e-10'
303  call exit(201)
304  endif
305  do j=1,3
306  xnoref(j)=xnoref(j)/dd
307  enddo
308  endif
309 !
310 ! updating the pointers iponor
311 !
312  nexp=nexp+1
313  do j=1,nel
314  if(ial(j).eq.1) then
315  ial(j)=2
316  if(ifix.eq.0) then
317  iponor(ipkonfa(iel(j))+jl(j))=ixfree
318  elseif(j.ne.jact) then
319  iponor(ipkonfa(iel(j))+jl(j))=
320  & iponor(ipkonfa(iel(jact))+jl(jact))
321  endif
322  endif
323  enddo
324 !
325 ! storing the normal in xnor and generating 3 nodes
326 ! for knor
327 !
328  if(ifix.eq.0) then
329  do j=1,3
330  xnor(ixfree+j)=xnoref(j)
331  enddo
332  ixfree=ixfree+3
333  endif
334 !
335  enddo
336  endif
337 !
338 ! write equations in file "jobname.equ"
339 !
340 ! check if node is a designvariable
341 !
342  if(nodedesiinv(i).eq.0) then
343 !
344 ! write equations in case nexp is greater or equal 3
345 !
346  if(nexp.ge.3) then
347  do j=1,3
348  write(20,106) 1
349  write(20,105) i,j,1
350  enddo
351 !
352 ! write equations in case nexp is 1
353 !
354  elseif(nexp.eq.1) then
355  j=1
356  do l=1,3
357  xnorloc(4-l,j)=xnor(ixfree+1-l)
358  sort(4-l)=dabs(xnor(ixfree+1-l))
359  nsort(4-l)=4-l
360  enddo
361  call dsort(sort,nsort,3,2)
362  write(20,106) 3
363  write(20,104) i,nsort(3),xnorloc(nsort(3),1),
364  & i,nsort(2),xnorloc(nsort(2),1),
365  & i,nsort(1),xnorloc(nsort(1),1)
366 !
367 ! write equations in case nexp is 2
368 !
369  elseif(nexp.eq.2) then
370  do j=1,nexp
371  k=j*3-3
372  do l=1,3
373  xnorloc(4-l,j)=xnor(ixfree+1-l-k)
374  enddo
375  enddo
376  if(i.eq.5457) then
377  ndet(1)=1
378  endif
379  ndet(1)=1
380  ndet(2)=2
381  ndet(3)=3
382  det(1)=dabs(xnorloc(1,1)*xnorloc(2,2)-
383  & xnorloc(1,2)*xnorloc(2,1))
384  det(2)=dabs(xnorloc(1,1)*xnorloc(3,2)-
385  & xnorloc(1,2)*xnorloc(3,1))
386  det(3)=dabs(xnorloc(2,1)*xnorloc(3,2)-
387  & xnorloc(2,2)*xnorloc(3,1))
388  call dsort(det,ndet,3,2)
389 
390  if(ndet(3).eq.1) then
391 ! if((dabs(xnorloc(1,1)).ge.dabs(xnorloc(2,1))).and.
392 ! & (dabs(xnorloc(2,2)).gt.1.d-5)) then
393  if((dabs(xnorloc(1,1)).gt.1.d-5).and.
394  & (dabs(xnorloc(2,2)).gt.1.d-5)) then
395  write(20,106) 2
396  write(20,107) i,1,xnorloc(1,1),
397  & i,2,xnorloc(2,1)
398  write(20,106) 2
399  write(20,107) i,2,xnorloc(2,2),
400  & i,1,xnorloc(1,2)
401  else
402  write(20,106) 2
403  write(20,107) i,2,xnorloc(2,1),
404  & i,1,xnorloc(1,1)
405  write(20,106) 2
406  write(20,107) i,1,xnorloc(1,2),
407  & i,2,xnorloc(2,2)
408  endif
409  elseif(ndet(3).eq.2) then
410 ! if((dabs(xnorloc(1,1)).ge.dabs(xnorloc(3,1))).and.
411 ! & (dabs(xnorloc(3,2)).gt.1.d-5)) then
412  if((dabs(xnorloc(1,1)).gt.1.d-5).and.
413  & (dabs(xnorloc(3,2)).gt.1.d-5)) then
414  write(20,106) 2
415  write(20,107) i,1,xnorloc(1,1),
416  & i,3,xnorloc(3,1)
417  write(20,106) 2
418  write(20,107) i,3,xnorloc(3,2),
419  & i,1,xnorloc(1,2)
420  else
421  write(20,106) 2
422  write(20,107) i,3,xnorloc(3,1),
423  & i,1,xnorloc(1,1)
424  write(20,106) 2
425  write(20,107) i,1,xnorloc(1,2),
426  & i,3,xnorloc(3,2)
427  endif
428  elseif(ndet(3).eq.3) then
429 ! if((dabs(xnorloc(2,1)).ge.dabs(xnorloc(3,1))).and.
430 ! & (dabs(xnorloc(3,2)).gt.1.d-5)) then
431  if((dabs(xnorloc(2,1)).gt.1.d-5).and.
432  & (dabs(xnorloc(3,2)).gt.1.d-5)) then
433  write(20,106) 2
434  write(20,107) i,2,xnorloc(2,1),
435  & i,3,xnorloc(3,1)
436  write(20,106) 2
437  write(20,107) i,3,xnorloc(3,2),
438  & i,2,xnorloc(2,2)
439  else
440  write(20,106) 2
441  write(20,107) i,3,xnorloc(3,1),
442  & i,2,xnorloc(2,1)
443  write(20,106) 2
444  write(20,107) i,2,xnorloc(2,2),
445  & i,3,xnorloc(3,2)
446  endif
447  endif
448  endif
449  endif
450 !
451  104 format(3(i10,",",i1,",",e20.13,","))
452  105 format(1(i10,",",i1,",",i1,","))
453  106 format(i1)
454  107 format(2(i10,",",i1,",",e20.13,","))
455 !
456  enddo
457  close(20)
458 !
459  return
460 !
461  100 write(*,*) '*ERROR in openfile: could not open file ',
462  & fnequ(1:len+4)
463  call exit(201)
464  101 write(*,*) '*ERROR in openfile: could not delete file ',
465  & fnequ(1:len+4)
466  call exit(201)
467 !
norshell8
subroutine norshell8(xi, et, xl, xnor)
Definition: norshell8.f:20
norshell4
subroutine norshell4(xi, et, xl, xnor)
Definition: norshell4.f:20
norshell6
subroutine norshell6(xi, et, xl, xnor)
Definition: norshell6.f:20
dsort
subroutine dsort(dx, iy, n, kflag)
Definition: dsort.f:6
norshell3
subroutine norshell3(xi, et, xl, xnor)
Definition: norshell3.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)