diff --git a/example/nearest_point_1d.f90 b/example/nearest_point_1d.f90 index d1a236d09..962bce256 100644 --- a/example/nearest_point_1d.f90 +++ b/example/nearest_point_1d.f90 @@ -48,7 +48,8 @@ program nearest_point_1d ! nearest_Xt: Corresponding parametric coordinates of the nearest point (optional) ! id: id of the nearest point (optional) call shape%nearest_point([4.5_rk, 4.5_rk, 5.0_rk], nearest_Xg, nearest_Xt, id) - print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,a,1x,g0)','Nearest point on the curve:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id + print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,a,1x,g0)',& + 'Nearest point on the curve:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id !----------------------------------------------------------------------------- ! Nearest point on the curve (Optimization) @@ -60,7 +61,8 @@ program nearest_point_1d ! nearest_Xt: Corresponding parametric coordinates of the nearest point ! nearest_Xg: Coordinates of the nearest point on the curve (optional) call shape%nearest_point2([4.5_rk, 4.5_rk, 5.0_rk], 1.0e-11_rk, 30, nearest_Xt, nearest_Xg) - print '(a,1x,g0,2x,g0,a,2x,g0,2x,g0)', 'Nearest point on the curve:', nearest_Xg, ' with parametric coordinates:', nearest_Xt + print '(a,1x,g0,2x,g0,a,2x,g0,2x,g0)',& + 'Nearest point on the curve:', nearest_Xg, ' with parametric coordinates:', nearest_Xt !----------------------------------------------------------------------------- ! Finalizing diff --git a/example/nearest_point_2d.f90 b/example/nearest_point_2d.f90 index 974543093..5de270744 100644 --- a/example/nearest_point_2d.f90 +++ b/example/nearest_point_2d.f90 @@ -42,7 +42,8 @@ program nearest_point_2d ! nearest_Xt: Corresponding parametric coordinates of the nearest point (optional) ! id: id of the nearest point (optional) call shape%nearest_point([1.3_rk, 1.0_rk, 1.999999999_rk], nearest_Xg, nearest_Xt, id) - print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,a,1x,g0)','Nearest point on the surface:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id + print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,a,1x,g0)',& + 'Nearest point on the surface:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id !----------------------------------------------------------------------------- ! Nearest point on the surface (Optimization) @@ -54,7 +55,8 @@ program nearest_point_2d ! nearest_Xt: Corresponding parametric coordinates of the nearest point ! nearest_Xg: Coordinates of the nearest point on the surface (optional) call shape%nearest_point2([1.3_rk, 1.0_rk, 1.999999999_rk], 1.0e-11_rk, 30, nearest_Xt, nearest_Xg) - print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0)', 'Nearest point on the surface:', nearest_Xg, ' with parametric coordinates:', nearest_Xt + print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0)',& + 'Nearest point on the surface:', nearest_Xg, ' with parametric coordinates:', nearest_Xt !----------------------------------------------------------------------------- ! Finalizing diff --git a/example/nearest_point_3d.f90 b/example/nearest_point_3d.f90 index 34f216c8b..41682e922 100644 --- a/example/nearest_point_3d.f90 +++ b/example/nearest_point_3d.f90 @@ -27,7 +27,11 @@ program nearest_point_3d !> The weights of the control points (Wc) are optional. Wc = [1.0_rk, 1.1_rk, 1.11_rk, 1.0_rk, 0.5_rk, 0.5_rk, 1.2_rk, 1.0_rk] - call shape%set(knot1=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], knot2=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], knot3=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk], Xc=Xc, Wc=Wc) + call shape%set(& + knot1=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk],& + knot2=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk],& + knot3=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk],& + Xc=Xc, Wc=Wc) !----------------------------------------------------------------------------- ! Creating the NURBS volume @@ -45,7 +49,8 @@ program nearest_point_3d ! nearest_Xt: Corresponding parametric coordinates of the nearest point (optional) ! id: id of the nearest point (optional) call shape%nearest_point([1.5_rk, 3.5_rk, 1.1_rk], nearest_Xg, nearest_Xt, id) - print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,g0,2x,a,1x,g0)','Nearest point on the volume:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id + print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,g0,2x,a,1x,g0)',& + 'Nearest point on the volume:', nearest_Xg, ' with parametric coordinates:', nearest_Xt, ' and id:', id !----------------------------------------------------------------------------- ! Nearest point on the volume (Optimization) @@ -57,7 +62,8 @@ program nearest_point_3d ! nearest_Xt: Corresponding parametric coordinates of the nearest point ! nearest_Xg: Coordinates of the nearest point on the volume (optional) call shape%nearest_point2([1.5_rk, 3.5_rk, 1.1_rk], 1.0e-11_rk, 500, nearest_Xt, nearest_Xg) - print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,g0)', 'Nearest point on the volume:', nearest_Xg, ' with parametric coordinates:', nearest_Xt + print '(a,1x,g0,2x,g0,2x,g0,a,2x,g0,2x,g0,2x,g0)',& + 'Nearest point on the volume:', nearest_Xg, ' with parametric coordinates:', nearest_Xt !----------------------------------------------------------------------------- ! Finalizing diff --git a/src/forcad_nurbs_curve.f90 b/src/forcad_nurbs_curve.f90 index 20fbd7428..f59c4c8f0 100644 --- a/src/forcad_nurbs_curve.f90 +++ b/src/forcad_nurbs_curve.f90 @@ -1806,7 +1806,8 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest obj = norm2(Xg - point_Xg) + 0.001_rk ! add a small number to avoid division by zero grad = dot_product((Xg-point_Xg)/obj, matmul(dTgc,this%Xc)) - hess = dot_product(matmul(dTgc,this%Xc) - (Xg-point_Xg)/obj*grad, matmul(dTgc,this%Xc))/obj + dot_product((Xg-point_Xg)/obj, matmul(d2Tgc,this%Xc)) + hess = dot_product(matmul(dTgc,this%Xc) - (Xg-point_Xg)/obj*grad, matmul(dTgc,this%Xc))/obj& + + dot_product((Xg-point_Xg)/obj, matmul(d2Tgc,this%Xc)) ! debug print '(i3,1x,e20.10,1x,e20.10)', k, xk, abs(grad) diff --git a/src/forcad_nurbs_surface.f90 b/src/forcad_nurbs_surface.f90 index f0cdcd606..259f4e524 100644 --- a/src/forcad_nurbs_surface.f90 +++ b/src/forcad_nurbs_surface.f90 @@ -2431,10 +2431,18 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest grad(1) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,1),this%Xc)) grad(2) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,2),this%Xc)) - hess(1,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2) ,1),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(1) ) / obj**2 - hess(2,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)+1:2*this%nc(1)*this%nc(2),1),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(1) ) / obj**2 - hess(1,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2) ,2),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(2) ) / obj**2 - hess(2,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)+1:2*this%nc(1)*this%nc(2),2),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(2) ) / obj**2 + hess(1,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,1),this%Xc)) + & + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2) ,1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(1) ) / obj**2 + hess(2,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,2),this%Xc)) + & + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)+1:2*this%nc(1)*this%nc(2),1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(1) ) / obj**2 + hess(1,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,1),this%Xc)) + & + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2) ,2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(2) ) / obj**2 + hess(2,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,2),this%Xc)) + & + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)+1:2*this%nc(1)*this%nc(2),2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(2) ) / obj**2 ! debug print '(i3,1x,2e20.10,1x,e20.10)', k, xk, norm2(grad) @@ -2776,10 +2784,18 @@ impure subroutine compute_d2Tgc_nurbs_2d_vector(Xt, knot1, knot2, degree, nc, ng d2Bi(1:nc(1)*nc(2) ,2) = kron(dB2, dB1) d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = kron(d2B2, B1) - d2Tgc(i,1:nc(1)*nc(2) ,1) = (d2Bi(1:nc(1)*nc(2) ,1)*Wc - 2.0_rk*dTgc(i,:,1)*dot_product(dBi(:,1),Wc) - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2) ,1),Wc)) / dot_product(Bi,Wc) - d2Tgc(i,nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1)*Wc - dTgc(i,:,1)*dot_product(dBi(:,2),Wc) - dTgc(i,:,2)*dot_product(dBi(:,1),Wc) - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1),Wc)) / dot_product(Bi,Wc) - d2Tgc(i,1:nc(1)*nc(2) ,2) = (d2Bi(1:nc(1)*nc(2) ,2)*Wc - dTgc(i,:,1)*dot_product(dBi(:,2),Wc) - dTgc(i,:,2)*dot_product(dBi(:,1),Wc) - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2) ,2),Wc)) / dot_product(Bi,Wc) - d2Tgc(i,nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2)*Wc - 2.0_rk*dTgc(i,:,2)*dot_product(dBi(:,2),Wc) - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i,1:nc(1)*nc(2) ,1) = & + (d2Bi(1:nc(1)*nc(2) ,1)*Wc - 2.0_rk*dTgc(i,:,1)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i,nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = & + (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1)*Wc - dTgc(i,:,1)*dot_product(dBi(:,2),Wc) - dTgc(i,:,2)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i,1:nc(1)*nc(2) ,2) = & + (d2Bi(1:nc(1)*nc(2) ,2)*Wc - dTgc(i,:,1)*dot_product(dBi(:,2),Wc) - dTgc(i,:,2)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i,nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = & + (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2)*Wc - 2.0_rk*dTgc(i,:,2)*dot_product(dBi(:,2),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2),Wc)) / dot_product(Bi,Wc) end do end subroutine !=============================================================================== @@ -2829,10 +2845,18 @@ impure subroutine compute_d2Tgc_nurbs_2d_scalar(Xt, knot1, knot2, degree, nc, Wc d2Bi(1:nc(1)*nc(2) ,2) = kron(dB2, dB1) d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = kron(d2B2, B1) - d2Tgc(1:nc(1)*nc(2) ,1) = (d2Bi(1:nc(1)*nc(2) ,1)*Wc - 2.0_rk*dTgc(:,1)*dot_product(dBi(:,1),Wc) - Tgc*dot_product(d2Bi(1:nc(1)*nc(2) ,1),Wc)) / dot_product(Bi,Wc) - d2Tgc(nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1)*Wc - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) - Tgc*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1),Wc)) / dot_product(Bi,Wc) - d2Tgc(1:nc(1)*nc(2) ,2) = (d2Bi(1:nc(1)*nc(2) ,2)*Wc - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) - Tgc*dot_product(d2Bi(1:nc(1)*nc(2) ,2),Wc)) / dot_product(Bi,Wc) - d2Tgc(nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2)*Wc - 2.0_rk*dTgc(:,2)*dot_product(dBi(:,2),Wc) - Tgc*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2),Wc)) / dot_product(Bi,Wc) + d2Tgc(1:nc(1)*nc(2) ,1) = & + (d2Bi(1:nc(1)*nc(2) ,1)*Wc - 2.0_rk*dTgc(:,1)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)+1:2*nc(1)*nc(2),1) = & + (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1)*Wc - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),1),Wc)) / dot_product(Bi,Wc) + d2Tgc(1:nc(1)*nc(2) ,2) = & + (d2Bi(1:nc(1)*nc(2) ,2)*Wc - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)+1:2*nc(1)*nc(2),2) = & + (d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2)*Wc - 2.0_rk*dTgc(:,2)*dot_product(dBi(:,2),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)+1:2*nc(1)*nc(2),2),Wc)) / dot_product(Bi,Wc) end subroutine !=============================================================================== diff --git a/src/forcad_nurbs_volume.f90 b/src/forcad_nurbs_volume.f90 index 284eb8192..31bce847f 100644 --- a/src/forcad_nurbs_volume.f90 +++ b/src/forcad_nurbs_volume.f90 @@ -205,7 +205,8 @@ pure subroutine compute_dTgc_bspline_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, end interface interface compute_d2Tgc - pure subroutine compute_d2Tgc_nurbs_3d_vector(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Wc, f_d2Tgc, f_dTgc, f_Tgc) + pure subroutine compute_d2Tgc_nurbs_3d_vector(& + f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_Wc, f_d2Tgc, f_dTgc, f_Tgc) import :: rk real(rk), intent(in), contiguous :: f_Xt(:,:) real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) @@ -218,7 +219,8 @@ pure subroutine compute_d2Tgc_nurbs_3d_vector(f_Xt, f_knot1, f_knot2, f_knot3, f real(rk), allocatable, intent(out) :: f_Tgc(:,:) end subroutine - pure subroutine compute_d2Tgc_bspline_3d_vector(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_d2Tgc, f_dTgc, f_Tgc) + pure subroutine compute_d2Tgc_bspline_3d_vector(& + f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_ng, f_d2Tgc, f_dTgc, f_Tgc) import :: rk real(rk), intent(in), contiguous :: f_Xt(:,:) real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) @@ -2949,15 +2951,33 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest grad(2) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,2),this%Xc)) grad(3) = dot_product((Xg-point_Xg)/obj, matmul(dTgc(:,3),this%Xc)) - hess(1,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(1) ) / obj**2 - hess(2,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(1) ) / obj**2 - hess(3,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(1) ) / obj**2 - hess(1,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(2) ) / obj**2 - hess(2,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(2) ) / obj**2 - hess(3,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(2) ) / obj**2 - hess(1,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(3) ) / obj**2 - hess(2,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(3) ) / obj**2 - hess(3,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(3) ) / obj**2 + hess(1,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(1) ) / obj**2 + hess(2,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(1) ) / obj**2 + hess(3,1) = ( dot_product(matmul(dTgc(:,1),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,1),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(1) ) / obj**2 + hess(1,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(2) ) / obj**2 + hess(2,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(2) ) / obj**2 + hess(3,2) = ( dot_product(matmul(dTgc(:,2),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,2),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(2) ) / obj**2 + hess(1,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,1),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(1:this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,1), this%Xc))*grad(3) ) / obj**2 + hess(2,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,2),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(this%nc(1)*this%nc(2)*this%nc(3)+1:2*this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,2), this%Xc))*grad(3) ) / obj**2 + hess(3,3) = ( dot_product(matmul(dTgc(:,3),this%Xc), matmul(dTgc(:,3),this%Xc)) + dot_product((Xg-point_Xg), & + matmul(d2Tgc(2*this%nc(1)*this%nc(2)*this%nc(3)+1:3*this%nc(1)*this%nc(2)*this%nc(3) ,3),this%Xc)) ) /obj & + - ( dot_product(Xg-point_Xg, matmul(dTgc(:,3), this%Xc))*grad(3) ) / obj**2 ! debug print '(i3,1x,3e20.10,1x,e20.10)', k, xk, norm2(grad) @@ -3518,15 +3538,33 @@ impure subroutine compute_d2Tgc_nurbs_3d_vector(Xt, knot1, knot2, knot3, degree, d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,dB2),B1) d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = kron(kron(d2B3,B2),B1) - d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,1) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,1)*Wc - 2.0_rk*dTgc(i, :,1)*dot_product(dBi(:,1),Wc) - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) - d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1)*Wc - dTgc(i, :,1)*dot_product(dBi(:,2),Wc) - dTgc(i, :,2)*dot_product(dBi(:,1),Wc) - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) - d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1)*Wc - dTgc(i, :,1)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,1),Wc) - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) - d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,2) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,2)*Wc - dTgc(i, :,1)*dot_product(dBi(:,2),Wc) - dTgc(i, :,2)*dot_product(dBi(:,1),Wc) - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) - d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2)*Wc - 2.0_rk*dTgc(i, :,2)*dot_product(dBi(:,2),Wc) - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) - d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2)*Wc - dTgc(i, :,2)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,2),Wc) - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) - d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,3) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,3)*Wc - dTgc(i, :,1)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,1),Wc) - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) - d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3)*Wc - dTgc(i, :,2)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,2),Wc) - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) - d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3)*Wc - 2.0_rk*dTgc(i, :,3)*dot_product(dBi(:,3),Wc) - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,1) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,1)*Wc & + - 2.0_rk*dTgc(i, :,1)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1)*Wc & + - dTgc(i, :,1)*dot_product(dBi(:,2),Wc) - dTgc(i, :,2)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1)*Wc & + - dTgc(i, :,1)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,2) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,2)*Wc & + - dTgc(i, :,1)*dot_product(dBi(:,2),Wc) - dTgc(i, :,2)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2)*Wc & + - 2.0_rk*dTgc(i, :,2)*dot_product(dBi(:,2),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2)*Wc & + - dTgc(i, :,2)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,2),Wc) & + - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 1:nc(1)*nc(2)*nc(3) ,3) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,3)*Wc & + - dTgc(i, :,1)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,1),Wc) & + - Tgc(i,:)*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3)*Wc & + - dTgc(i, :,2)*dot_product(dBi(:,3),Wc) - dTgc(i, :,3)*dot_product(dBi(:,2),Wc) & + - Tgc(i,:)*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(i, 2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3)*Wc & + - 2.0_rk*dTgc(i, :,3)*dot_product(dBi(:,3),Wc) & + - Tgc(i,:)*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) end do end subroutine !=============================================================================== @@ -3586,15 +3624,33 @@ impure subroutine compute_d2Tgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = kron(kron(dB3,dB2),B1) d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = kron(kron(d2B3,B2),B1) - d2Tgc(1:nc(1)*nc(2)*nc(3) ,1) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,1)*Wc - 2.0_rk*dTgc(:,1)*dot_product(dBi(:,1),Wc) - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) - d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1)*Wc - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) - d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1)*Wc - dTgc(:,1)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,1),Wc) - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) - d2Tgc(1:nc(1)*nc(2)*nc(3) ,2) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,2)*Wc - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) - d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2)*Wc - 2.0_rk*dTgc(:,2)*dot_product(dBi(:,2),Wc) - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) - d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2)*Wc - dTgc(:,2)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,2),Wc) - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) - d2Tgc(1:nc(1)*nc(2)*nc(3) ,3) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,3)*Wc - dTgc(:,1)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,1),Wc) - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) - d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3)*Wc - dTgc(:,2)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,2),Wc) - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) - d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3)*Wc - 2.0_rk*dTgc(:,3)*dot_product(dBi(:,3),Wc) - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(1:nc(1)*nc(2)*nc(3) ,1) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,1)*Wc & + - 2.0_rk*dTgc(:,1)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1)*Wc & + - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1)*Wc & + - dTgc(:,1)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,1),Wc)) / dot_product(Bi,Wc) + d2Tgc(1:nc(1)*nc(2)*nc(3) ,2) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,2)*Wc & + - dTgc(:,1)*dot_product(dBi(:,2),Wc) - dTgc(:,2)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2)*Wc & + - 2.0_rk*dTgc(:,2)*dot_product(dBi(:,2),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2)*Wc & + - dTgc(:,2)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,2),Wc) & + - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,2),Wc)) / dot_product(Bi,Wc) + d2Tgc(1:nc(1)*nc(2)*nc(3) ,3) = (d2Bi(1:nc(1)*nc(2)*nc(3) ,3)*Wc & + - dTgc(:,1)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,1),Wc) & + - Tgc*dot_product(d2Bi(1:nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3)*Wc & + - dTgc(:,2)*dot_product(dBi(:,3),Wc) - dTgc(:,3)*dot_product(dBi(:,2),Wc) & + - Tgc*dot_product(d2Bi(nc(1)*nc(2)*nc(3)+1:2*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) + d2Tgc(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3) = (d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3)*Wc & + - 2.0_rk*dTgc(:,3)*dot_product(dBi(:,3),Wc) & + - Tgc*dot_product(d2Bi(2*nc(1)*nc(2)*nc(3)+1:3*nc(1)*nc(2)*nc(3) ,3),Wc)) / dot_product(Bi,Wc) end subroutine !=============================================================================== diff --git a/test/fdm_curve.f90 b/test/fdm_curve.f90 index da429e683..325d8b0a8 100644 --- a/test/fdm_curve.f90 +++ b/test/fdm_curve.f90 @@ -8,7 +8,8 @@ program fdm_test_curve real(rk), allocatable :: Xc(:,:), Wc(:) !! Arrays for control points and weights real(rk) :: knot(6) !! Array for knot vector real(rk) :: Xtp, tol, Xt, Xtm - real(rk), allocatable :: Tgc(:), dTgc(:), Tgcp(:), dTgcp(:), Tgcm(:), dTgcm(:), CFD(:), BFD(:), FFD(:), d2Tgc(:), d2Tgcp(:), d2Tgcm(:) + real(rk), allocatable :: Tgc(:), dTgc(:), Tgcp(:), dTgcp(:), Tgcm(:), dTgcm(:), & + CFD(:), BFD(:), FFD(:), d2Tgc(:), d2Tgcp(:), d2Tgcm(:) !----------------------------------------------------------------------------- ! Setting up the NURBS curve diff --git a/test/fdm_surface.f90 b/test/fdm_surface.f90 index ed7f850e8..593328bf1 100644 --- a/test/fdm_surface.f90 +++ b/test/fdm_surface.f90 @@ -14,7 +14,8 @@ program fdm_test_surface !----------------------------------------------------------------------------- ! Setting up the NURBS surface !----------------------------------------------------------------------------- - Wc = [1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk] + Wc = [1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk,& + 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk, 1.0_rk, 0.5_rk, 1.0_rk, 0.2_rk] call surface%set_tetragon(L=[5.0_rk, 8.0_rk], nc=[4,4], Wc=Wc) !-----------------------------------------------------------------------------