proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
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
12c ----------------------------------------------------------------------
13c... 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
25c ------------------------------------------------------------------
26c... Element number,hess_flag
27c... 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
31c... Vector of Local basis function vales at (u_hat, v_hat), local and
32c 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
41c... Local Variables
42c
43c 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
47c 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
52c NURBS coordinates, counters for loops
53 integer ni, nj, nk, i, j, k, icount, aa
54
55c temporary variables
56 real*8 tmp
57
58c ----------------------------------------------------------------------
59
60
61c 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
70c 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
76c 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
84c 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
97c... 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
102c... 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
120c 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
138c
139c... Now calculate gradients.
140c
141
142c 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
187c
188c.... comPte the inverse of deformation gradient
189c
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
250c... 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
280c... 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
314c... global Hessian
315
316c... 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
343c... RHS of the matrix eQation for the second derivatives of bases.
344c... 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
365c... 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
409ccc write(*,*)"lhs"
410ccc write(*,*) locLHS(:,:)
411ccc write(*,*)"rhs"
412ccc write(*,*) shhessl(1,:)
413ccc write(*,*)"ans"
414
415c... (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
434ccc write(*,*) shhessl(1,:)
435ccc write(*,*) "end"
436
437c... 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
472c --------------VARIABLE DECLARATIONS-------------------------------
473c... 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
475c 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
480c ------------------------------------------------------------------
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
539c 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
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
w
#define w(x)
Definition jf.h:22
Generated on for proteus by doxygen 1.18.0