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

Go to the source code of this file.

Functions/Subroutines

subroutine chksurf (lakon, kon, ipkon, neigh, ipneigh, co, itypflag, node, icont, iscount, angmax)
 

Function/Subroutine Documentation

◆ chksurf()

subroutine chksurf ( character*8, dimension(*)  lakon,
integer, dimension(*)  kon,
integer, dimension(*)  ipkon,
integer, dimension(2,*)  neigh,
integer, dimension(*)  ipneigh,
real*8, dimension(3,*)  co,
integer  itypflag,
integer  node,
integer  icont,
integer  iscount,
real*8  angmax 
)
21 !
22 ! icont=1: element surfaces adjacent to a surface node have normal
23 ! vectors which have an angle of less than 10 degree
24 ! -> free surface assumed
25 ! icont=0: -> edge assumed
26 !
27 ! also counts the free surfaces adjacent to a node
28 !
29 ! author: Sascha Merz
30 !
31  implicit none
32 !
33  integer kon(*),ipkon(*),ielem,i,j,k,indexe,
34  & neigh(2,*),ipneigh(*),index,m,nvertex,itypflag,isurf,node,index1,
35  & ielem1,ncount,ntos8h(3,8),ntos4tet(3,4),iston8h(4,6),
36  & isnode, isidx,ifreesur(3),iston20h(8,6),iston10tet(6,4),
37  & iscount,lnod,icont,iston4tet(3,4)
38 !
39  real*8 co(3,*),angle,shpder8q(2,4,8),shpder6tri(2,3,6),
40  & vectors(3,3),vlen(2),lastvec(3),angtmp,xl(3,8),
41  & shpder4q(2,4,4),shpder3tri(2,3,3),angmax
42 !
43 ! ntosX(j,k) returns the three surface id's j for the corner node k
44 ! for the element surfaces adjacent to the node
45 !
46  data ntos8h /1,3,6,1,3,4,1,4,5,1,5,6,2,3,6,2,3,4,2,4,5,2,5,6/
47 !
48  data ntos4tet /1,2,4,1,2,3,1,3,4,2,3,4/
49 !
50 ! istonX(j,k) returns the nodes j of the element surface k
51 !
52  data iston8h /1,2,3,4,5,8,7,6,1,5,6,2,2,6,7,3,3,7,8,4,4,8,5,1/
53 !
54  data iston20h /1,2,3,4,9,10,11,12,5,8,7,6,16,15,14,13,
55  & 1,5,6,2,17,13,18,9,2,6,7,3,18,14,19,10,
56  & 3,7,8,4,19,15,20,11,4,8,5,1,20,16,17,12/
57 !
58  data iston4tet /1,2,3,1,4,2,2,4,3,3,4,1/
59 !
60  data iston10tet /1,2,3,5,6 ,7,1,4,2,8 ,9,5,
61  & 2,4,3,9,10,6,3,4,1,10,8,7/
62 !
63 ! shpder8q contains the first derivative of the shape functions
64 ! of a 8 node quadrilateral element with shpder8q(i,j,k) where i
65 ! can be 1 for the xi-derivative or 2 for the eta-derivative j
66 ! can be 1-4 for the location in the corner nodes to be evaluated
67 ! for the element nodes k
68 !
69  data shpder8q /
70  & -1.5d0,-1.5d0, 0.5d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.5d0,
71  & -0.5d0, 0.0d0, 1.5d0,-1.5d0, 0.0d0, 0.5d0, 0.0d0, 0.0d0,
72  & 0.0d0, 0.0d0, 0.0d0,-0.5d0, 1.5d0, 1.5d0,-0.5d0, 0.0d0,
73  & 0.0d0,-0.5d0, 0.0d0, 0.0d0, 0.5d0, 0.0d0,-1.5d0, 1.5d0,
74  & 2.0d0, 0.0d0,-2.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,
75  & 0.0d0, 0.0d0, 0.0d0, 2.0d0, 0.0d0,-2.0d0, 0.0d0, 0.0d0,
76  & 0.0d0, 0.0d0, 0.0d0, 0.0d0,-2.0d0, 0.0d0, 2.0d0, 0.0d0,
77  & 0.0d0, 2.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,-2.0d0/
78 !
79 ! same as above for a 4 node linear quadrilateral element
80 !
81  data shpder4q /
82  & -0.5d0,-0.5d0,-0.5d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0,-0.5d0,
83  & 0.5d0, 0.0d0, 0.5d0,-0.5d0, 0.0d0,-0.5d0, 0.0d0, 0.0d0,
84  & 0.0d0, 0.0d0, 0.0d0, 0.5d0, 0.5d0, 0.5d0, 0.5d0, 0.0d0,
85  & 0.0d0, 0.5d0, 0.0d0, 0.0d0,-0.5d0, 0.0d0,-0.5d0, 0.5d0/
86 !
87 ! same as above for a 6 node quadratic triangular element
88 !
89  data shpder6tri /
90  & -3.0d0,-3.0d0, 1.0d0, 1.0d0, 1.0d0, 1.0d0,
91  & -1.0d0, 0.0d0, 3.0d0, 0.0d0,-1.0d0, 0.0d0,
92  & 0.0d0,-1.0d0, 0.0d0,-1.0d0, 0.0d0, 3.0d0,
93  & 4.0d0, 0.0d0,-4.0d0,-4.0d0, 0.0d0, 0.0d0,
94  & 0.0d0, 0.0d0, 0.0d0, 4.0d0, 4.0d0, 0.0d0,
95  & 0.0d0, 4.0d0, 0.0d0, 0.0d0,-4.0d0,-4.0d0/
96 !
97 ! same as above for a 3 node linear triangular element
98 !
99  data shpder3tri /
100  & -1.0d0,-1.0d0,-1.0d0,-1.0d0,-1.0d0,-1.0d0,
101  & 1.0d0, 0.0d0, 1.0d0, 0.0d0, 1.0d0, 0.0d0,
102  & 0.0d0, 1.0d0, 0.0d0, 1.0d0, 0.0d0, 1.0d0/
103 !
104  character*8 lakon(*)
105 !
106  index=ipneigh(node)
107  icont=1
108  iscount=0
109  angmax=0.d0
110 !
111  do
112  if(index.eq.0) exit
113  ielem=neigh(1,index)
114 !
115  if(lakon(ielem)(1:5).eq.'C3D20'.and.itypflag.eq.1) then
116  nvertex=8
117  elseif(lakon(ielem)(1:5).eq.'C3D10'.and.itypflag.eq.2) then
118  nvertex=4
119  elseif(lakon(ielem)(1:4).eq.'C3D8'.and.itypflag.eq.3) then
120  nvertex=8
121  elseif(lakon(ielem)(1:4).eq.'C3D4'.and.itypflag.eq.4) then
122  nvertex=4
123  else
124  index=neigh(2,index)
125  cycle
126  endif
127 !
128 ! find the index of the node in the element
129 !
130  indexe=ipkon(ielem)
131  do m=1,nvertex
132  if(kon(indexe+m).eq.node) exit
133  enddo
134 !
135 ! the local node number is m
136 !
137 ! now every surface has to be checked
138 !
139  do j=1,3
140  ifreesur(j)=0
141  enddo
142 !
143  do isurf=1,3
144 !
145 ! finding the global node numbers of the
146 ! nodes of the surface
147 !
148  if(nvertex.eq.4) then
149 !
150 ! isidx: index of the surface neighbouring the node
151 !
152  isidx=ntos4tet(isurf,m)
153  elseif(nvertex.eq.8) then
154  isidx=ntos8h(isurf,m)
155  endif
156 !
157 ! find out, if there is any element neighbouring 'node',
158 ! which has also those nodes (-> surface is within volume)
159 !
160  index1=ipneigh(node)
161  do
162  if(index1.eq.0) exit
163  ielem1=neigh(1,index1)
164  if(
165  & .not.(
166  & lakon(ielem1)(1:5).eq.'C3D20'.and.itypflag.eq.1
167  & .or.
168  & lakon(ielem1)(1:5).eq.'C3D10'.and.itypflag.eq.2
169  & .or.
170  & lakon(ielem1)(1:4).eq.'C3D8'.and.itypflag.eq.3
171  & .or.
172  & lakon(ielem1)(1:4).eq.'C3D4'.and.itypflag.eq.4
173  & )
174  & .or.ielem.eq.ielem1
175  & ) then
176  index1=neigh(2,index1)
177  cycle
178  endif
179 !
180 ! check every corner node in the element
181 !
182  ncount=0
183  do k=1,3
184  if(nvertex.eq.4) then
185  isnode=kon(indexe+iston4tet(k,isidx))
186  elseif(nvertex.eq.8) then
187  isnode=kon(indexe+iston8h(k,isidx))
188  endif
189  do j=1,nvertex
190  if(kon(ipkon(ielem1)+j).eq.isnode)
191  & ncount=ncount+1
192  enddo
193  enddo
194 !
195  if(ncount.eq.3) then
196 !
197 ! surface isurf is not a free surface
198 !
199  ifreesur(isurf)=1
200  endif
201 !
202  index1=neigh(2,index1)
203  enddo
204  enddo
205 !
206  do isurf=1,3
207  if(ifreesur(isurf).eq.0) then
208  iscount=iscount+1
209  do i=1,3
210  do j=1,3
211  vectors(j,i)=0.d0
212  enddo
213  enddo
214 !
215 ! free surface: find out local node number
216 ! of the 'surface element' neighbouring the
217 ! node to be evaluated
218 !
219  if(nvertex.eq.8) then
220  isidx=ntos8h(isurf,m)
221  do j=1,4
222  if( (isidx.eq.1.and.m.eq.1)
223  & .or.(isidx.eq.2.and.m.eq.5)
224  & .or.(isidx.eq.3.and.m.eq.1)
225  & .or.(isidx.eq.4.and.m.eq.2)
226  & .or.(isidx.eq.5.and.m.eq.3)
227  & .or.(isidx.eq.6.and.m.eq.4)) then
228  lnod=1
229  elseif( (isidx.eq.1.and.m.eq.2)
230  & .or.(isidx.eq.2.and.m.eq.8)
231  & .or.(isidx.eq.3.and.m.eq.5)
232  & .or.(isidx.eq.4.and.m.eq.6)
233  & .or.(isidx.eq.5.and.m.eq.7)
234  & .or.(isidx.eq.6.and.m.eq.8)) then
235  lnod=2
236  elseif( (isidx.eq.1.and.m.eq.3)
237  & .or.(isidx.eq.2.and.m.eq.7)
238  & .or.(isidx.eq.3.and.m.eq.6)
239  & .or.(isidx.eq.4.and.m.eq.7)
240  & .or.(isidx.eq.5.and.m.eq.8)
241  & .or.(isidx.eq.6.and.m.eq.5)) then
242  lnod=3
243  elseif( (isidx.eq.1.and.m.eq.4)
244  & .or.(isidx.eq.2.and.m.eq.6)
245  & .or.(isidx.eq.3.and.m.eq.2)
246  & .or.(isidx.eq.4.and.m.eq.3)
247  & .or.(isidx.eq.5.and.m.eq.4)
248  & .or.(isidx.eq.6.and.m.eq.1)) then
249  lnod=4
250  endif
251  enddo
252 !
253 c do k=1,8
254 c write(*,*) 'node',node,' nodes',
255 c & kon(indexe+iston20h(k,isidx)),'lnod',lnod
256 c enddo
257 !
258  if(itypflag.eq.1) then
259 !
260 ! get coordinates of the 2d-element nodes
261 !
262  do k=1,8
263  do j=1,3
264  xl(j,k)=co(j,kon(indexe+iston20h(k,isidx)))
265  enddo
266  enddo
267 !
268 ! vectors(j,i) (i=1,2) is the j-derivative for the
269 ! coordinates i.
270 !
271  do k=1,8
272  do i=1,3
273  do j=1,2
274  vectors(j,i)=vectors(j,i)+
275  & xl(i,k)*shpder8q(j,lnod,k)
276  enddo
277  enddo
278  enddo
279 !
280  elseif(itypflag.eq.3) then
281  do k=1,4
282  do j=1,3
283  xl(j,k)=co(j,kon(indexe+iston8h(k,isidx)))
284  enddo
285  enddo
286  do k=1,4
287  do i=1,3
288  do j=1,2
289  vectors(j,i)=vectors(j,i)+
290  & xl(i,k)*shpder4q(j,lnod,k)
291  enddo
292  enddo
293  enddo
294  endif
295 !
296  elseif(nvertex.eq.4) then
297  isidx=ntos4tet(isurf,m)
298  do j=1,3
299  if( (isidx.eq.1.and.m.eq.1)
300  & .or.(isidx.eq.2.and.m.eq.1)
301  & .or.(isidx.eq.3.and.m.eq.2)
302  & .or.(isidx.eq.4.and.m.eq.3)) then
303  lnod=1
304  elseif( (isidx.eq.1.and.m.eq.2)
305  & .or.(isidx.eq.2.and.m.eq.4)
306  & .or.(isidx.eq.3.and.m.eq.4)
307  & .or.(isidx.eq.4.and.m.eq.4)) then
308  lnod=2
309  elseif( (isidx.eq.1.and.m.eq.3)
310  & .or.(isidx.eq.2.and.m.eq.2)
311  & .or.(isidx.eq.3.and.m.eq.3)
312  & .or.(isidx.eq.4.and.m.eq.1)) then
313  lnod=3
314  endif
315  enddo
316 !
317  if(itypflag.eq.2) then
318  do k=1,6
319  do j=1,3
320  xl(j,k)=co(j,kon(indexe+iston10tet(k,isidx)))
321  enddo
322  enddo
323  do k=1,6
324  do i=1,3
325  do j=1,2
326  vectors(j,i)=vectors(j,i)+
327  & xl(i,k)*shpder6tri(j,lnod,k)
328  enddo
329  enddo
330  enddo
331  elseif(itypflag.eq.4) then
332  do k=1,3
333  do j=1,3
334  xl(j,k)=co(j,kon(indexe+iston4tet(k,isidx)))
335  enddo
336  enddo
337  do k=1,3
338  do i=1,3
339  do j=1,2
340  vectors(j,i)=vectors(j,i)+
341  & xl(i,k)*shpder3tri(j,lnod,k)
342  enddo
343  enddo
344  enddo
345  endif
346 !
347  endif
348 !
349 ! vectors(3,i) is the normal vector of the surface in the
350 ! evaluated node 'node'
351 !
352  vectors(3,1)=vectors(1,2)*vectors(2,3)
353  & -vectors(1,3)*vectors(2,2)
354  vectors(3,2)=vectors(1,3)*vectors(2,1)
355  & -vectors(1,1)*vectors(2,3)
356  vectors(3,3)=vectors(1,1)*vectors(2,2)
357  & -vectors(1,2)*vectors(2,1)
358  vlen(2)=dsqrt(vectors(3,1)*vectors(3,1)
359  & +vectors(3,2)*vectors(3,2)
360  & +vectors(3,3)*vectors(3,3))
361 !
362  if(iscount.gt.1) then
363  angtmp=dabs((vectors(3,1)*lastvec(1)
364  & +vectors(3,2)*lastvec(2)
365  & +vectors(3,3)*lastvec(3))
366  & /(vlen(1)*vlen(2)))
367  if(angtmp.lt.1.d0) then
368  angle=57.29577951d0*dacos(angtmp)
369  if(angle.gt.angmax) angmax=angle
370  else
371  angle=0.d0
372  endif
373 !
374 ! if the angle between the normal vectors of two
375 ! surfaces is greater than 10 degree, than it is
376 ! assumed that an edge (dicontinuity) is present
377 !
378  if(angle.ge.10.d0) then
379  icont=0
380  endif
381  endif
382 !
383  do j=1,3
384  lastvec(j)=vectors(3,j)
385  enddo
386  vlen(1)=vlen(2)
387 !
388  endif
389  enddo
390  index=neigh(2,index)
391  enddo
392  return
Hosted by OpenAircraft.com, (Michigan UAV, LLC)