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

Go to the source code of this file.

Functions/Subroutines

subroutine basis (x, y, z, xo, yo, zo, nx, ny, nz, planfa, ifatet, nktet, netet, field, nfield, cotet, kontyp, ipkon, kon, iparent, xp, yp, zp, value, ratio, iselect, nselect, istartset, iendset, ialset, imastset, ielemnr, nterms, konl)
 

Function/Subroutine Documentation

◆ basis()

subroutine basis ( real*8, dimension(*), intent(in)  x,
real*8, dimension(*), intent(in)  y,
real*8, dimension(*), intent(in)  z,
real*8, dimension(*), intent(in)  xo,
real*8, dimension(*), intent(in)  yo,
real*8, dimension(*), intent(in)  zo,
integer, dimension(*), intent(in)  nx,
integer, dimension(*), intent(in)  ny,
integer, dimension(*), intent(in)  nz,
real*8, dimension(4,*), intent(in)  planfa,
integer, dimension(4,*), intent(in)  ifatet,
integer, intent(in)  nktet,
integer, intent(in)  netet,
real*8, dimension(nfield,nktet), intent(in)  field,
integer, intent(in)  nfield,
real*8, dimension(3,*), intent(in)  cotet,
integer, dimension(*), intent(in)  kontyp,
integer, dimension(*), intent(in)  ipkon,
integer, dimension(*), intent(in)  kon,
integer, dimension(*), intent(in)  iparent,
real*8, intent(in)  xp,
real*8, intent(in)  yp,
real*8, intent(in)  zp,
real*8, dimension(nfield), intent(inout)  value,
real*8, dimension(20), intent(inout)  ratio,
integer, dimension(nselect), intent(in)  iselect,
integer, intent(in)  nselect,
integer, dimension(*), intent(in)  istartset,
integer, dimension(*), intent(in)  iendset,
integer, dimension(*), intent(in)  ialset,
integer, intent(in)  imastset,
integer, dimension(*), intent(in)  ielemnr,
integer, intent(inout)  nterms,
integer, dimension(20), intent(inout)  konl 
)
25 !
26  implicit none
27 !
28 ! 100 nearest nodes: node(100),idummy1(100),idummy2(100),iparentel(100)
29 !
30  integer ifatet(4,*),id,node(100),near,nx(*),ny(*),nz(*),
31  & ifs,iface,i,konl(20),nfield,nktet,ielement,ielmax,netet,k,
32  & j,iselect(nselect),nselect,kontyp(*),ipkon(*),kon(*),iparent(*),
33  & nterms,indexe,nelem,konl_opt(20),idummy1(100),idummy2(100),
34  & iparentel(100),nparentel,kflag,ii,inside,nterms_opt,
35  & istartset(*),iendset(*),ialset(*),imastset,ielemnr(*),
36  & ielementnr,nlength,nearset
37 !
38  real*8 cotet(3,*),planfa(4,*),dface,dfacemax,tolerance,
39  & dist,field(nfield,nktet),ratio(20),pneigh(3,20),pnode(3),xi,et,
40  & ze,dist_opt,ratio_opt(20),xp,yp,zp,x(*),y(*),z(*),xo(*),yo(*),
41  & zo(*),value(nfield)
42 !
43  intent(in) x,y,z,xo,yo,zo,nx,ny,nz,
44  & planfa,ifatet,nktet,netet,field,nfield,
45  & cotet,kontyp,ipkon,kon,iparent,
46  & xp,yp,zp,iselect,nselect,
47  & istartset,iendset,ialset,imastset,ielemnr
48 !
49  intent(inout) konl,value,nterms,ratio
50 !
51  tolerance=1.d-6
52 !
53 ! look for the global element encompassing point (xp,yp,zp)
54 !
55  do ii=1,3
56 !
57 ! three-level algorithm
58 ! - start with the search for the nearest neighboring element
59 ! - if this element is not the right one search for the 10
60 ! nearest neighbors
61 ! - if no fitting element was found, search for the 100 nearest
62 ! neighbors
63 !
64  if(ii.eq.1) then
65  near=1
66  elseif(ii.eq.2) then
67  near=10
68  else
69  near=100
70  endif
71  call near3d(xo,yo,zo,x,y,z,nx,ny,nz,xp,yp,zp,netet,
72  & node,near)
73 !
74  inside=0
75  nearset=0
76  ielmax=0
77 !
78  do i=1,near
79  ielement=node(i)
80 c write(*,*) 'basis tetrahedron ',ielement,iparent(ielement)
81 !
82 ! check whether the element belongs to the right set,
83 ! if a set is specified
84 !
85 ! ielement is a tetrahedral element number
86 ! iparent(ielement) is a running number from 1 to the
87 ! total number of elements
88 ! ielemnr(iparent(ielement)) is the real element number
89 ! (there can be gaps in the numbering)
90 !
91  if(imastset.ne.0) then
92  ielementnr=ielemnr(iparent(ielement))
93  nlength=iendset(imastset)-istartset(imastset)+1
94  call nident(ialset(istartset(imastset)),ielementnr,
95  & nlength,id)
96  if(id.le.0) cycle
97  if(ialset(istartset(imastset)+id-1).ne.ielementnr) cycle
98  nearset=nearset+1
99  node(nearset)=node(i)
100  else
101  nearset=near
102  endif
103 
104  dface=0.d0
105  do j=1,4
106 !
107 ! the face number in ifatet is negative, if the equation
108 ! of the face plane is such, that its value in the
109 ! remaining node of the tetrahedron is negative
110 !
111  ifs=ifatet(j,ielement)
112  iface=abs(ifs)
113  dist=planfa(1,iface)*xp+planfa(2,iface)*yp
114  & +planfa(3,iface)*zp+planfa(4,iface)
115  if(dist*ifs.lt.-1.d-10*iface) then
116  dface=dface+dist*ifs/iface
117  endif
118  enddo
119 c write(*,*) 'basis dface ',dface
120  if(dface.gt.-1.d-10) then
121  inside=1
122  exit
123  endif
124 c if(i.eq.1) then
125  if(ielmax.eq.0) then
126  dfacemax=dface
127  ielmax=ielement
128  else
129  if(dface.gt.dfacemax) then
130  ielmax=ielement
131  dfacemax=dface
132  endif
133  endif
134  enddo
135 !
136 ! if a global element set was defined, check whether an
137 ! appropriate element was found
138 !
139  if(imastset.ne.0) then
140  if(nearset.eq.0) then
141  if(ii.lt.3) then
142  cycle
143  else
144  write(*,*) '*ERROR: no suitable global element found'
145  write(*,*) ' for location (',xp,yp,zp,')'
146  call exit(201)
147  endif
148  endif
149  endif
150 !
151 ! if no element was found, the element with the smallest
152 ! dfacemax (in absolute value; summed distance) is taken
153 !
154  if(inside.eq.0) then
155  ielement=ielmax
156  endif
157 !
158  nelem=iparent(ielement)
159 c write(*,*) 'basis element ',nelem
160  if(kontyp(nelem).eq.1) then
161  nterms=8
162  elseif(kontyp(nelem).eq.2) then
163  nterms=6
164  elseif(kontyp(nelem).eq.3) then
165  nterms=4
166  elseif(kontyp(nelem).eq.4) then
167  nterms=20
168  elseif(kontyp(nelem).eq.5) then
169  nterms=15
170  elseif(kontyp(nelem).eq.6) then
171  nterms=10
172  else
173  nterms=0
174  endif
175  indexe=ipkon(nelem)
176 !
177 ! modify order of connectivity ccx <-> cgx
178 !
179  if(kontyp(nelem).eq.4) then
180  do i=1,12
181  konl(i)=kon(indexe+i)
182  enddo
183  do i=13,16
184  konl(i+4)=kon(indexe+i)
185  enddo
186  do i=17,20
187  konl(i-4)=kon(indexe+i)
188  enddo
189  elseif(kontyp(nelem).eq.5) then
190  do i=1,9
191  konl(i)=kon(indexe+i)
192  enddo
193  do i=10,12
194  konl(i+3)=kon(indexe+i)
195  enddo
196  do i=13,15
197  konl(i-3)=kon(indexe+i)
198  enddo
199  else
200  do i=1,nterms
201  konl(i)=kon(indexe+i)
202  enddo
203  endif
204 !
205 ! nodes of master element
206 !
207  do i=1,nterms
208  do j=1,3
209  pneigh(j,i)=cotet(j,konl(i))
210  enddo
211  enddo
212 !
213 ! slave node
214 !
215  pnode(1)=xp
216  pnode(2)=yp
217  pnode(3)=zp
218 !
219 ! attaching slave node to master element
220 !
221  if(nterms.ne.0) then
222  call attach_3d(pneigh,pnode,nterms,ratio,dist,xi,et,ze)
223  endif
224 !
225 ! checking the parent elements of the "best" tetrahedra
226 ! in case the distance between slave node and location of
227 ! interpolation exceeds "tolerance"
228 !
229  if(dist.gt.tolerance) then
230  do i=1,nterms
231  konl_opt(i)=konl(i)
232  ratio_opt(i)=ratio(i)
233  enddo
234  nterms_opt=nterms
235  dist_opt=dist
236 !
237 ! sorting the parent elements
238 !
239  do i=1,nearset
240  idummy1(i)=iparent(node(i))
241  enddo
242  kflag=1
243  call isortii(idummy1,idummy2,nearset,kflag)
244  nparentel=0
245  do i=1,nearset
246  if(idummy1(i).eq.nelem) cycle
247  if(i.eq.1) then
248  iparentel(1)=idummy1(1)
249  nparentel=1
250  else
251  if(idummy1(i).ne.idummy1(i-1)) then
252  nparentel=nparentel+1
253  iparentel(nparentel)=idummy1(i)
254  endif
255  endif
256  enddo
257 !
258  do k=1,nparentel
259 !
260 ! slave node
261 !
262  pnode(1)=xp
263  pnode(2)=yp
264  pnode(3)=zp
265  nelem=iparentel(k)
266  if(kontyp(nelem).eq.1) then
267  nterms=8
268  elseif(kontyp(nelem).eq.2) then
269  nterms=6
270  elseif(kontyp(nelem).eq.3) then
271  nterms=4
272  elseif(kontyp(nelem).eq.4) then
273  nterms=20
274  elseif(kontyp(nelem).eq.5) then
275  nterms=15
276  elseif(kontyp(nelem).eq.6) then
277  nterms=10
278  else
279  nterms=0
280  endif
281  indexe=ipkon(nelem)
282 !
283 ! modify order of connectivity ccx <-> cgx
284 !
285  if(kontyp(nelem).eq.4) then
286  do i=1,12
287  konl(i)=kon(indexe+i)
288  enddo
289  do i=13,16
290  konl(i+4)=kon(indexe+i)
291  enddo
292  do i=17,20
293  konl(i-4)=kon(indexe+i)
294  enddo
295  elseif(kontyp(nelem).eq.5) then
296  do i=1,9
297  konl(i)=kon(indexe+i)
298  enddo
299  do i=10,12
300  konl(i+3)=kon(indexe+i)
301  enddo
302  do i=13,15
303  konl(i-3)=kon(indexe+i)
304  enddo
305  else
306  do i=1,nterms
307  konl(i)=kon(indexe+i)
308  enddo
309  endif
310 !
311 ! nodes of master element
312 !
313  do i=1,nterms
314  do j=1,3
315  pneigh(j,i)=cotet(j,konl(i))
316  enddo
317  enddo
318 !
319 ! attaching slave node to master element
320 !
321  if(nterms.ne.0) then
322  call attach_3d(pneigh,pnode,nterms,ratio,dist,
323  & xi,et,ze)
324  endif
325 !
326 ! check whether the present element yields better results
327 !
328  if(dist.lt.dist_opt) then
329  do i=1,nterms
330  konl_opt(i)=konl(i)
331  ratio_opt(i)=ratio(i)
332  enddo
333  nterms_opt=nterms
334  dist_opt=dist
335  endif
336  if(dist.lt.tolerance) exit
337  enddo
338 !
339 ! storing the optimal configuration
340 !
341  nterms=nterms_opt
342  do i=1,nterms
343  konl(i)=konl_opt(i)
344  ratio(i)=ratio_opt(i)
345  enddo
346  endif
347 !
348  if((ii.eq.3).or.(dist.lt.tolerance)) exit
349 !
350  enddo
351 !
352 ! interpolating the fields
353 !
354  do k=1,nselect
355  i=iselect(k)
356  value(k)=0.d0
357  do j=1,nterms
358  value(k)=value(k)+ratio(j)*field(i,konl(j))
359 c write(*,*) 'basis j',j,konl(j),field(i,konl(j))
360  enddo
361 c write(*,*) 'basis select',k,i,value(k)
362  enddo
363 !
364  return
near3d
subroutine near3d(xo, yo, zo, x, y, z, nx, ny, nz, xp, yp, zp, n, neighbor, k)
Definition: near3d.f:20
attach_3d
subroutine attach_3d(pneigh, pnode, nterms, ratio, dist, xil, etl, zel)
Definition: attach_3d.f:20
isortii
subroutine isortii(ix, iy, n, kflag)
Definition: isortii.f:6
dist
static double * dist
Definition: radflowload.c:42
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
Hosted by OpenAircraft.com, (Michigan UAV, LLC)