proteus  1.8.1
C/C++/Fortran libraries
NURBSshapes.f
Go to the documentation of this file.
1 
2 
3  subroutine eval_shape_3d(
4  & NSD,P,Q,R,NSHL,MCP,NCP,OCP,NEL,
5  & U_KNOT,V_KNOT,W_KNOT,B_NET,
6  & ni,nj,nk,u_hat,v_hat,w_hat,nders,
7  & shl,shgradl,shgradg,shhessg,dxidx,DetJ
8  & )
9 
10  implicit none
11 
12 c ----------------------------------------------------------------------
13 c... Define stuff that would normally be in common/adjkeep
14  integer NSD
15 
16  integer P,Q,R,NSHL
17  integer MCP,NCP,OCP,NEL
18 
19  real*8 u_knot(mcp+p+1),v_knot(ncp+q+1),w_knot(ocp+r+1)
20  real*8 b_net(mcp,ncp,ocp,nsd+1)
21  !real*8 B_NET(P+1,Q+1,R+1,NCP,OCP,NSD+1)
22 
23  real DetJ
24 
25 c ------------------------------------------------------------------
26 c... Element number,hess_flag
27 c... u and v coordinates of integration point in parent element
28  real*8 u_hat, v_hat, w_hat, du, dv, dw
29  integer e,nders
30 
31 c... Vector of Local basis function vales at (u_hat, v_hat), local and
32 c global gradients.
33 
34 
35  real*8 shl(nshl),shgradl(nshl,nsd),shgradg(nshl,nsd),
36  & shhessg(nshl,nsd,nsd), tempshl(nshl),tempshgradl(nshl,nsd),
37  & shhessl(nshl,6), tempshhessl(nshl,6)
38 
39  real*8 dxdxi(nsd,nsd), dxidx(nsd,nsd),dxdxixj(nsd,6), loclhs(6,6)
40 
41 c... Local Variables
42 c
43 c 1D nonrational basis functions and derivs in u and v
44 
45  real*8 n(nders+1,p+1), m(nders+1,q+1),o(nders+1,r+1)
46 
47 c u and v coordinates of integration point, denominator and derivative sums
48  real*8 u, v, w, denom_sum, derv_sum_u, derv_sum_v,
49  & derv_sum_w, derv_sum_uu,
50  & derv_sum_uv, derv_sum_uw,derv_sum_vv, derv_sum_vw,derv_sum_ww
51 
52 c NURBS coordinates, counters for loops
53  integer ni, nj, nk, i, j, k, icount, aa
54 
55 c temporary variables
56  real*8 tmp
57 
58 c ----------------------------------------------------------------------
59 
60 
61 c Get u and v coordinates of integration point
62 
63  u = ((u_knot(ni+1)-u_knot(ni))*u_hat +
64  & u_knot(ni+1) + u_knot(ni))/2d+0
65  v = ((v_knot(nj+1)-v_knot(nj))*v_hat +
66  & v_knot(nj+1) + v_knot(nj))/2d+0
67  w = ((w_knot(nk+1)-w_knot(nk))*w_hat +
68  & w_knot(nk+1) + w_knot(nk))/2d+0
69 
70 c Get knot span sizes
71 
72  du = u_knot(ni+1)-u_knot(ni)
73  dv = v_knot(nj+1)-v_knot(nj)
74  dw = w_knot(nk+1)-w_knot(nk)
75 
76 c Evalate 1D shape functions and derivatives each direction
77 
78 
79  call dersbasisfuns(ni,p,mcp,u,nders,u_knot,n) ! calculate in u dir.
80  call dersbasisfuns(nj,q,ncp,v,nders,v_knot,m) ! calculate in v dir.
81  call dersbasisfuns(nk,r,ocp,w,nders,w_knot,o) ! calculate in v dir.
82 
83 
84 c Form basis functions and derivatives dR/du and dR/dv
85 
86  icount = 0
87  denom_sum = 0d+0
88  derv_sum_u = 0d+0
89  derv_sum_v = 0d+0
90  derv_sum_w = 0d+0
91 
92  do k = 0, r
93  do j = 0, q
94  do i = 0, p
95  icount = icount+1
96 
97 c... basis functions
98  shl(icount) = n(1,p+1-i)*m(1,q+1-j)*o(1,r+1-k)*
99  & b_net(ni-i,nj-j,nk-k,nsd+1)
100  denom_sum = denom_sum + shl(icount)
101 
102 c... derivatives
103  shgradl(icount,1) =
104  & n(2,p+1-i)*m(1,q+1-j)*o(1,r+1-k)*
105  & b_net(ni-i,nj-j,nk-k,nsd+1) ! u
106  derv_sum_u = derv_sum_u + shgradl(icount,1)
107  shgradl(icount,2) =
108  & n(1,p+1-i)*m(2,q+1-j)*o(1,r+1-k)*
109  & b_net(ni-i,nj-j,nk-k,nsd+1) ! v
110  derv_sum_v = derv_sum_v + shgradl(icount,2)
111  shgradl(icount,3) =
112  & n(1,p+1-i)*m(1,q+1-j)*o(2,r+1-k)*
113  & b_net(ni-i,nj-j,nk-k,nsd+1) ! w
114  derv_sum_w = derv_sum_w + shgradl(icount,3)
115 
116  enddo
117  enddo
118  enddo
119 
120 c Divide through by denominator
121 
122  tempshl = shl
123  tempshgradl = shgradl
124 
125 
126  do i = 1,nshl
127  shgradl(i,1) = shgradl(i,1)/denom_sum -
128  & (shl(i)*derv_sum_u)/(denom_sum**2)
129  shgradl(i,2) = shgradl(i,2)/denom_sum -
130  & (shl(i)*derv_sum_v)/(denom_sum**2)
131  shgradl(i,3) = shgradl(i,3)/denom_sum -
132  & (shl(i)*derv_sum_w)/(denom_sum**2)
133  shl(i) = shl(i)/denom_sum
134  enddo
135 
136 
137 
138 c
139 c... Now calculate gradients.
140 c
141 
142 c calculate dx/dxi
143 
144 
145 
146  dxdxi = 0d+0
147  icount = 0
148 
149  do k = 0, r
150  do j = 0, q
151  do i = 0, p
152  icount = icount + 1
153 
154  dxdxi(1,1) = dxdxi(1,1) +
155  & b_net(ni-i,nj-j,nk-k,1) *
156  & shgradl(icount,1)
157  dxdxi(1,2) = dxdxi(1,2) +
158  & b_net(ni-i,nj-j,nk-k,1) *
159  & shgradl(icount,2)
160  dxdxi(1,3) = dxdxi(1,3) +
161  & b_net(ni-i,nj-j,nk-k,1) *
162  & shgradl(icount,3)
163  dxdxi(2,1) = dxdxi(2,1) +
164  & b_net(ni-i,nj-j,nk-k,2) *
165  & shgradl(icount,1)
166  dxdxi(2,2) = dxdxi(2,2) +
167  & b_net(ni-i,nj-j,nk-k,2) *
168  & shgradl(icount,2)
169  dxdxi(2,3) = dxdxi(2,3) +
170  & b_net(ni-i,nj-j,nk-k,2) *
171  & shgradl(icount,3)
172  dxdxi(3,1) = dxdxi(3,1) +
173  & b_net(ni-i,nj-j,nk-k,3) *
174  & shgradl(icount,1)
175  dxdxi(3,2) = dxdxi(3,2) +
176  & b_net(ni-i,nj-j,nk-k,3) *
177  & shgradl(icount,2)
178  dxdxi(3,3) = dxdxi(3,3) +
179  & b_net(ni-i,nj-j,nk-k,3) *
180  & shgradl(icount,3)
181 
182  enddo
183  enddo
184  enddo
185 
186 
187 c
188 c.... comPte the inverse of deformation gradient
189 c
190 
191  dxidx(1,1) = dxdxi(2,2) * dxdxi(3,3)
192  & - dxdxi(3,2) * dxdxi(2,3)
193  dxidx(1,2) = dxdxi(3,2) * dxdxi(1,3)
194  & - dxdxi(1,2) * dxdxi(3,3)
195  dxidx(1,3) = dxdxi(1,2) * dxdxi(2,3)
196  & - dxdxi(1,3) * dxdxi(2,2)
197  tmp = 1d+0 / ( dxidx(1,1) * dxdxi(1,1)
198  & + dxidx(1,2) * dxdxi(2,1)
199  & + dxidx(1,3) * dxdxi(3,1) )
200  dxidx(1,1) = dxidx(1,1) * tmp
201  dxidx(1,2) = dxidx(1,2) * tmp
202  dxidx(1,3) = dxidx(1,3) * tmp
203  dxidx(2,1) = (dxdxi(2,3) * dxdxi(3,1)
204  & - dxdxi(2,1) * dxdxi(3,3)) * tmp
205  dxidx(2,2) = (dxdxi(1,1) * dxdxi(3,3)
206  & - dxdxi(3,1) * dxdxi(1,3)) * tmp
207  dxidx(2,3) = (dxdxi(2,1) * dxdxi(1,3)
208  & - dxdxi(1,1) * dxdxi(2,3)) * tmp
209  dxidx(3,1) = (dxdxi(2,1) * dxdxi(3,2)
210  & - dxdxi(2,2) * dxdxi(3,1)) * tmp
211  dxidx(3,2) = (dxdxi(3,1) * dxdxi(1,2)
212  & - dxdxi(1,1) * dxdxi(3,2)) * tmp
213  dxidx(3,3) = (dxdxi(1,1) * dxdxi(2,2)
214  & - dxdxi(1,2) * dxdxi(2,1)) * tmp
215 
216  detj = 1d+0/tmp ! Note that DetJ resides in common
217 
218 
219 
220 
221  do i = 1, nshl
222 
223  shgradg(i,1) = shgradl(i,1) * dxidx(1,1) +
224  & shgradl(i,2) * dxidx(2,1) +
225  & shgradl(i,3) * dxidx(3,1)
226  shgradg(i,2) = shgradl(i,1) * dxidx(1,2) +
227  & shgradl(i,2) * dxidx(2,2) +
228  & shgradl(i,3) * dxidx(3,2)
229  shgradg(i,3) = shgradl(i,1) * dxidx(1,3) +
230  & shgradl(i,2) * dxidx(2,3) +
231  & shgradl(i,3) * dxidx(3,3)
232 
233  enddo
234 
235 
236  if (nders.eq.2) then
237 
238  icount = 0
239  derv_sum_uu = 0d+0
240  derv_sum_uv = 0d+0
241  derv_sum_uw = 0d+0
242  derv_sum_vv = 0d+0
243  derv_sum_vw = 0d+0
244  derv_sum_ww = 0d+0
245  do k = 0, r
246  do j = 0,q
247  do i = 0,p
248  icount = icount+1
249 
250 c... 2nd derivatives
251  tempshhessl(icount,1) =
252  & n(3,p+1-i)*m(1,q+1-j)*o(1,r+1-k)*
253  & b_net(ni-i,nj-j,nk-k,nsd+1) ! ,uu
254  derv_sum_uu = derv_sum_uu + tempshhessl(icount,1)
255  tempshhessl(icount,2) =
256  & n(2,p+1-i)*m(2,q+1-j)*o(1,r+1-k)*
257  & b_net(ni-i,nj-j,nk-k,nsd+1) ! ,uv
258  derv_sum_uv = derv_sum_uv + tempshhessl(icount,2)
259  tempshhessl(icount,3) =
260  & n(2,p+1-i)*m(1,q+1-j)*o(2,r+1-k)*
261  & b_net(ni-i,nj-j,nk-k,nsd+1) ! ,uw
262  derv_sum_uw = derv_sum_uw + tempshhessl(icount,3)
263  tempshhessl(icount,4) =
264  & n(1,p+1-i)*m(3,q+1-j)*o(1,r+1-k)*
265  & b_net(ni-i,nj-j,nk-k,nsd+1) ! ,vv
266  derv_sum_vv = derv_sum_vv + tempshhessl(icount,4)
267  tempshhessl(icount,5) =
268  & n(1,p+1-i)*m(2,q+1-j)*o(2,r+1-k)*
269  & b_net(ni-i,nj-j,nk-k,nsd+1) ! ,vw
270  derv_sum_vw = derv_sum_vw + tempshhessl(icount,5)
271  tempshhessl(icount,6) =
272  & n(1,p+1-i)*m(1,q+1-j)*o(3,r+1-k)*
273  & b_net(ni-i,nj-j,nk-k,nsd+1) ! ,ww
274  derv_sum_ww = derv_sum_ww + tempshhessl(icount,6)
275  enddo
276  enddo
277  enddo
278 
279 
280 c... local Hessian
281  do i = 1,nshl
282  shhessl(i,1) = tempshhessl(i,1)/denom_sum -
283  & tempshgradl(i,1)*derv_sum_u/(denom_sum**2) -
284  & ((tempshgradl(i,1)*derv_sum_u+tempshl(i)*derv_sum_uu)/
285  & (denom_sum**2) - 2d+0*tempshl(i)*derv_sum_u*derv_sum_u/
286  & (denom_sum**3))
287  shhessl(i,2) = tempshhessl(i,2)/denom_sum -
288  & tempshgradl(i,1)*derv_sum_v/(denom_sum**2) -
289  & ((tempshgradl(i,2)*derv_sum_u+tempshl(i)*derv_sum_uv)/
290  & (denom_sum**2) - 2d+0*tempshl(i)*derv_sum_u*derv_sum_v/
291  & (denom_sum**3))
292  shhessl(i,3) = tempshhessl(i,3)/denom_sum -
293  & tempshgradl(i,1)*derv_sum_w/(denom_sum**2) -
294  & ((tempshgradl(i,3)*derv_sum_u+tempshl(i)*derv_sum_uw)/
295  & (denom_sum**2) - 2d+0*tempshl(i)*derv_sum_u*derv_sum_w/
296  & (denom_sum**3))
297  shhessl(i,4) = tempshhessl(i,4)/denom_sum -
298  & tempshgradl(i,2)*derv_sum_v/(denom_sum**2) -
299  & ((tempshgradl(i,2)*derv_sum_v+tempshl(i)*derv_sum_vv)/
300  & (denom_sum**2) - 2d+0*tempshl(i)*derv_sum_v*derv_sum_v/
301  & (denom_sum**3))
302  shhessl(i,5) = tempshhessl(i,5)/denom_sum -
303  & tempshgradl(i,2)*derv_sum_w/(denom_sum**2) -
304  & ((tempshgradl(i,3)*derv_sum_v+tempshl(i)*derv_sum_vw)/
305  & (denom_sum**2) - 2d+0*tempshl(i)*derv_sum_v*derv_sum_w/
306  & (denom_sum**3))
307  shhessl(i,6) = tempshhessl(i,6)/denom_sum -
308  & tempshgradl(i,3)*derv_sum_w/(denom_sum**2) -
309  & ((tempshgradl(i,3)*derv_sum_w+tempshl(i)*derv_sum_ww)/
310  & (denom_sum**2) - 2d+0*tempshl(i)*derv_sum_w*derv_sum_w/
311  & (denom_sum**3))
312  enddo
313 
314 c... global Hessian
315 
316 c... Second derivatives of the geometrical map
317 
318  dxdxixj = 0d+0
319  icount = 0
320 
321  do k = 0, r
322  do j = 0, q
323  do i = 0, p
324  icount = icount + 1
325 
326  dxdxixj(1,:) = dxdxixj(1,:)+
327  & b_net(ni-i,nj-j,nk-k,1)*
328  & shhessl(icount,:)
329 
330  dxdxixj(2,:) = dxdxixj(2,:)+
331  & b_net(ni-i,nj-j,nk-k,2)*
332  & shhessl(icount,:)
333 
334  dxdxixj(3,:) = dxdxixj(3,:)+
335  & b_net(ni-i,nj-j,nk-k,3)*
336  & shhessl(icount,:)
337 
338 
339  enddo
340  enddo
341  enddo
342 
343 c... RHS of the matrix eQation for the second derivatives of bases.
344 c... Reuse local hess. array
345 
346  shhessl(:,1) = shhessl(:,1) - shgradg(:,1)*dxdxixj(1,1) -
347  & shgradg(:,2)*dxdxixj(2,1) - shgradg(:,3)*dxdxixj(3,1)
348 
349  shhessl(:,2) = shhessl(:,2) - shgradg(:,1)*dxdxixj(1,2) -
350  & shgradg(:,2)*dxdxixj(2,2) - shgradg(:,3)*dxdxixj(3,2)
351 
352  shhessl(:,3) = shhessl(:,3) - shgradg(:,1)*dxdxixj(1,3) -
353  & shgradg(:,2)*dxdxixj(2,3) - shgradg(:,3)*dxdxixj(3,3)
354 
355  shhessl(:,4) = shhessl(:,4) - shgradg(:,1)*dxdxixj(1,4) -
356  & shgradg(:,2)*dxdxixj(2,4) - shgradg(:,3)*dxdxixj(3,4)
357 
358  shhessl(:,5) = shhessl(:,5) - shgradg(:,1)*dxdxixj(1,5) -
359  & shgradg(:,2)*dxdxixj(2,5) - shgradg(:,3)*dxdxixj(3,5)
360 
361  shhessl(:,6) = shhessl(:,6) - shgradg(:,1)*dxdxixj(1,6) -
362  & shgradg(:,2)*dxdxixj(2,6) - shgradg(:,3)*dxdxixj(3,6)
363 
364 
365 c... LHS (6x6, same for every basis function)
366 
367  loclhs(1,1) = dxdxi(1,1)*dxdxi(1,1)
368  loclhs(1,2) = 2d+0*dxdxi(1,1)*dxdxi(2,1)
369  loclhs(1,3) = 2d+0*dxdxi(1,1)*dxdxi(3,1)
370  loclhs(1,4) = dxdxi(2,1)*dxdxi(2,1)
371  loclhs(1,5) = 2d+0*dxdxi(2,1)*dxdxi(3,1)
372  loclhs(1,6) = dxdxi(3,1)*dxdxi(3,1)
373 
374  loclhs(2,1) = dxdxi(1,1)*dxdxi(1,2)
375  loclhs(2,2) = dxdxi(1,1)*dxdxi(2,2) + dxdxi(1,2)*dxdxi(2,1)
376  loclhs(2,3) = dxdxi(1,1)*dxdxi(3,2) + dxdxi(1,2)*dxdxi(3,1)
377  loclhs(2,4) = dxdxi(2,1)*dxdxi(2,2)
378  loclhs(2,5) = dxdxi(2,1)*dxdxi(3,2) + dxdxi(2,2)*dxdxi(3,1)
379  loclhs(2,6) = dxdxi(3,1)*dxdxi(3,2)
380 
381  loclhs(3,1) = dxdxi(1,1)*dxdxi(1,3)
382  loclhs(3,2) = dxdxi(1,1)*dxdxi(2,3) + dxdxi(1,3)*dxdxi(2,1)
383  loclhs(3,3) = dxdxi(1,1)*dxdxi(3,3) + dxdxi(1,3)*dxdxi(3,1)
384  loclhs(3,4) = dxdxi(2,1)*dxdxi(2,3)
385  loclhs(3,5) = dxdxi(2,1)*dxdxi(3,3) + dxdxi(2,3)*dxdxi(3,1)
386  loclhs(3,6) = dxdxi(3,1)*dxdxi(3,3)
387 
388  loclhs(4,1) = dxdxi(1,2)*dxdxi(1,2)
389  loclhs(4,2) = 2d+0*dxdxi(1,2)*dxdxi(2,2)
390  loclhs(4,3) = 2d+0*dxdxi(1,2)*dxdxi(3,2)
391  loclhs(4,4) = dxdxi(2,2)*dxdxi(2,2)
392  loclhs(4,5) = 2d+0*dxdxi(2,2)*dxdxi(3,2)
393  loclhs(4,6) = dxdxi(3,2)*dxdxi(3,2)
394 
395  loclhs(5,1) = dxdxi(1,2)*dxdxi(1,3)
396  loclhs(5,2) = dxdxi(1,2)*dxdxi(2,3) + dxdxi(1,3)*dxdxi(2,2)
397  loclhs(5,3) = dxdxi(1,2)*dxdxi(3,3) + dxdxi(1,3)*dxdxi(3,2)
398  loclhs(5,4) = dxdxi(2,2)*dxdxi(2,3)
399  loclhs(5,5) = dxdxi(2,2)*dxdxi(3,3) + dxdxi(2,3)*dxdxi(3,2)
400  loclhs(5,6) = dxdxi(3,2)*dxdxi(3,3)
401 
402  loclhs(6,1) = dxdxi(1,3)*dxdxi(1,3)
403  loclhs(6,2) = 2d+0*dxdxi(1,3)*dxdxi(2,3)
404  loclhs(6,3) = 2d+0*dxdxi(1,3)*dxdxi(3,3)
405  loclhs(6,4) = dxdxi(2,3)*dxdxi(2,3)
406  loclhs(6,5) = 2d+0*dxdxi(2,3)*dxdxi(3,3)
407  loclhs(6,6) = dxdxi(3,3)*dxdxi(3,3)
408 
409 ccc write(*,*)"lhs"
410 ccc write(*,*) locLHS(:,:)
411 ccc write(*,*)"rhs"
412 ccc write(*,*) shhessl(1,:)
413 ccc write(*,*)"ans"
414 
415 c... (6x6) - Gaussian elimination
416 
417  do k = 1,6
418  do i = k+1,6
419  tmp = loclhs(i,k)/loclhs(k,k)
420  do j = k+1,6
421  loclhs(i,j) = loclhs(i,j) - tmp*loclhs(k,j)
422  enddo
423  shhessl(:,i) = shhessl(:,i) - tmp*shhessl(:,k)
424  enddo
425  enddo
426 
427  do i = 6,1,-1
428  do j = i+1,6
429  shhessl(:,i) = shhessl(:,i) - loclhs(i,j)*shhessl(:,j)
430  enddo
431  shhessl(:,i) = shhessl(:,i)/loclhs(i,i)
432  enddo
433 
434 ccc write(*,*) shhessl(1,:)
435 ccc write(*,*) "end"
436 
437 c... Assign to global hessian of basis functions
438 
439  shhessg(:,1,1) = shhessl(:,1)
440  shhessg(:,1,2) = shhessl(:,2)
441  shhessg(:,1,3) = shhessl(:,3)
442 
443  shhessg(:,2,1) = shhessl(:,2)
444  shhessg(:,2,2) = shhessl(:,4)
445  shhessg(:,2,3) = shhessl(:,5)
446 
447  shhessg(:,3,1) = shhessl(:,3)
448  shhessg(:,3,2) = shhessl(:,5)
449  shhessg(:,3,3) = shhessl(:,6)
450 
451  endif
452 
453  dxidx(1,:) = dxidx(1,:)*2d+0/du
454  dxidx(2,:) = dxidx(2,:)*2d+0/dv
455  dxidx(3,:) = dxidx(3,:)*2d+0/dw
456 
457 
458  detj = abs(detj)
459 
460  return
461  end
462 
463 
464 
465 
466 
467 !#######################################################################
468  subroutine dersbasisfuns(i,pl,ml,u,nders,u_knotl,ders)
469 
470  implicit none
471 
472 c --------------VARIABLE DECLARATIONS-------------------------------
473 c... knot span, degree of curve, number of control points, counters
474  integer i, pl, ml, j, r, k, j1, j2, s1, s2, rk, pk,nders
475 c parameter vale, vector of knots, derivative matrix
476  real*8 u, u_knotl(pl+ml+1), ders(nders+1,pl+1), ndu(pl+1,pl+1),d
477 
478  real*8 left(pl+1), right(pl+1), saved, temp, a(2,pl+1)
479 
480 c ------------------------------------------------------------------
481 
482  ndu(1,1) = 1
483  do j = 1,pl
484  left(j+1) = u - u_knotl(i+1-j)
485  right(j+1) = u_knotl(i+j) - u
486  saved = 0
487  do r = 0,j-1
488  ndu(j+1,r+1) = right(r+2) + left(j-r+1)
489  temp = ndu(r+1,j)/ndu(j+1,r+1)
490  ndu(r+1,j+1) = saved + right(r+2)*temp
491  saved = left(j-r+1)*temp
492  enddo
493  ndu(j+1,j+1) = saved
494  enddo
495 
496  ! load basis functions
497  do j = 0,pl
498  ders(1,j+1) = ndu(j+1,pl+1)
499  enddo
500  ! comPte derivatives
501  do r = 0,pl ! loop over function index
502  s1 = 0
503  s2 = 1 ! alternate rows in array a
504  a(1,1) = 1
505  ! loop to comPte kth derivative
506  do k = 1,nders
507  d = 0d+0
508  rk = r-k
509  pk = pl-k
510  if (r >= k) then
511  a(s2+1,1) = a(s1+1,1)/ndu(pk+2,rk+1)
512  d = a(s2+1,1)*ndu(rk+1,pk+1)
513  endif
514  if (rk >= -1) then
515  j1 = 1
516  else
517  j1 = -rk
518  endif
519  if ((r-1) <= pk) then
520  j2 = k-1
521  else
522  j2 = pl-r
523  endif
524  do j = j1,j2
525  a(s2+1,j+1) = (a(s1+1,j+1) - a(s1+1,j))/ndu(pk+2,rk+j+1)
526  d = d + a(s2+1,j+1)*ndu(rk+j+1,pk+1)
527  enddo
528  if (r <= pk) then
529  a(s2+1,k+1) = -a(s1+1,k)/ndu(pk+2,r+1)
530  d = d + a(s2+1,k+1)*ndu(r+1,pk+1)
531  endif
532  ders(k+1,r+1) = d
533  j = s1
534  s1 = s2
535  s2 = j ! switch rows
536  enddo
537  enddo
538 
539 c Multiply through by the correct factors
540  r = pl
541  do k = 1,nders
542  do j = 0,pl
543  ders(k+1,j+1) = ders(k+1,j+1)*r
544  enddo
545  r = r*(pl-k)
546  enddo
547 
548  return
549  end
dersbasisfuns
subroutine dersbasisfuns(i, pl, ml, u, nders, u_knotl, ders)
Definition: NURBSshapes.f:469
w
#define w(x)
Definition: jf.h:22
eval_shape_3d
subroutine eval_shape_3d(NSD, P, Q, R, NSHL, MCP, NCP, OCP, NEL, U_KNOT, V_KNOT, W_KNOT, B_NET, ni, nj, nk, u_hat, v_hat, w_hat, nders, shl, shgradl, shgradg, shhessg, dxidx, DetJ)
Definition: NURBSshapes.f:9
Generated on Fri Jul 1 2022 11:17:39 for proteus by doxygen 1.8.20