diff --git a/CHANGELOG.md b/CHANGELOG.md index 0b0418b0c..42fd904a2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,20 @@ All notable changes to [ForCAD](https://github.com/gha3mi/forcad) will be docume The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), +## [Unreleased] + +### Added + +- Added `det`, `inv`, `dyad` and `gauss_leg` procedures in the `forcad_utils` module. +- Added procedure `set1a` to the `nurbs_curve` derived type. +- Added procedure `set4` to the `nurbs_curve` NURBS derived type. +- Added optional input variable `elem` to `derivative_scalar` procedures. +- Added `ansatz` procedures to compute shape functions, derivatives of shape functions and (dV, dA, dL). +- Added `cmp_length()` to compute the length of a NURBS curve. +- Added `cmp_area()` to compute the area of a NURBS surface. +- Added `cmp_volume()` to compute the volume of a NURBS volume. +- Added examples for `cmp_length()`, `cmp_area()`, and `cmp_volume()`. + ## [0.5.1] ### Changed diff --git a/example/cmp_area.f90 b/example/cmp_area.f90 new file mode 100644 index 000000000..daca832bd --- /dev/null +++ b/example/cmp_area.f90 @@ -0,0 +1,23 @@ +program compute_area + + use forcad + + implicit none + + type(nurbs_surface) :: shape + real(rk) :: area + real(rk) :: Xc(4,3) + + Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk] + Xc(2,:) = [2.0_rk, 0.0_rk, 0.0_rk] + Xc(3,:) = [0.0_rk, 2.0_rk, 0.0_rk] + Xc(4,:) = [2.0_rk, 2.0_rk, 0.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],& + Xc=Xc) + + call shape%cmp_area(area) + print*, area +end program diff --git a/example/cmp_length.f90 b/example/cmp_length.f90 new file mode 100644 index 000000000..42e89adaa --- /dev/null +++ b/example/cmp_length.f90 @@ -0,0 +1,20 @@ +program compute_length + + use forcad + + implicit none + + type(nurbs_curve) :: shape + real(rk) :: length + real(rk) :: Xc(2,3) + + Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk] + Xc(2,:) = [2.0_rk, 0.0_rk, 0.0_rk] + + call shape%set(& + knot=[0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk],& + Xc=Xc) + + call shape%cmp_length(length) + print*, length +end program diff --git a/example/cmp_volume.f90 b/example/cmp_volume.f90 new file mode 100644 index 000000000..60efc1449 --- /dev/null +++ b/example/cmp_volume.f90 @@ -0,0 +1,28 @@ +program compute_volume + + use forcad + + implicit none + + type(nurbs_volume) :: shape + real(rk) :: volume + real(rk) :: Xc(8,3) + + Xc(1,:) = [0.0_rk, 0.0_rk, 0.0_rk] + Xc(2,:) = [2.0_rk, 0.0_rk, 0.0_rk] + Xc(3,:) = [0.0_rk, 2.0_rk, 0.0_rk] + Xc(4,:) = [2.0_rk, 2.0_rk, 0.0_rk] + Xc(5,:) = [0.0_rk, 0.0_rk, 2.0_rk] + Xc(6,:) = [2.0_rk, 0.0_rk, 2.0_rk] + Xc(7,:) = [0.0_rk, 2.0_rk, 2.0_rk] + Xc(8,:) = [2.0_rk, 2.0_rk, 2.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) + + call shape%cmp_volume(volume) + print*, volume +end program diff --git a/fpm.toml b/fpm.toml index 38855122f..4f32ed8cc 100644 --- a/fpm.toml +++ b/fpm.toml @@ -26,6 +26,7 @@ forimage = {git="https://github.com/gha3mi/forimage.git"} forcolormap = {git="https://github.com/vmagnin/forcolormap.git"} fortime = {git="https://github.com/gha3mi/fortime.git"} openmp = "*" +stdlib = "*" [dev-dependencies] forunittest = {git="https://github.com/gha3mi/forunittest"} @@ -158,4 +159,19 @@ main = "nearest_point_3d.f90" [[example]] name = "nearest_point_2d_bench" source-dir = "example" -main = "nearest_point_2d_bench.f90" \ No newline at end of file +main = "nearest_point_2d_bench.f90" + +[[example]] +name = "cmp_length" +source-dir = "example" +main = "cmp_length.f90" + +[[example]] +name = "cmp_area" +source-dir = "example" +main = "cmp_area.f90" + +[[example]] +name = "cmp_volume" +source-dir = "example" +main = "cmp_volume.f90" \ No newline at end of file diff --git a/src/forcad_nurbs_curve.f90 b/src/forcad_nurbs_curve.f90 index 6153670e8..34744fcf1 100644 --- a/src/forcad_nurbs_curve.f90 +++ b/src/forcad_nurbs_curve.f90 @@ -5,7 +5,7 @@ module forcad_nurbs_curve use forcad_utils, only: rk, basis_bspline, elemConn_C0, compute_multiplicity, compute_knot_vector, basis_bspline_der,& insert_knot_A_5_1, findspan, elevate_degree_A_5_9, remove_knots_A_5_8, & - elemConn_Cn, unique, rotation + elemConn_Cn, unique, rotation, dyad, gauss_leg implicit none @@ -29,9 +29,11 @@ module forcad_nurbs_curve integer, allocatable, private :: elemConn(:,:) !! IGA element connectivity contains procedure :: set1 !!> Set knot vector, control points and weights for the NURBS curve object + procedure :: set1a procedure :: set2 !!> Set NURBS curve using nodes of parameter space, degree, continuity, control points and weights procedure :: set3 !!> Set Bezier or Rational Bezier curve using control points and weights - generic :: set => set1, set2, set3 !!> Set NURBS curve + procedure :: set4 !!> Set NURBS curve using degree, number of control points, control points and weights + generic :: set => set1, set1a, set2, set3, set4 !!> Set NURBS curve procedure :: create !!> Generate geometry points procedure :: cmp_Xg !!> Compute geometry points procedure, private :: get_Xc_all !!> Get all control points @@ -90,6 +92,8 @@ module forcad_nurbs_curve procedure :: show !!> Show the NURBS object using PyVista procedure :: nearest_point !!> Find the nearest point on the NURBS curve (Approximation) procedure :: nearest_point2 !!> Find the nearest point on the NURBS curve (Minimization - Newton's method) + procedure :: ansatz !!> Compute the shape functions, derivative of shape functions and dL + procedure :: cmp_length !!> Compute the length of the NURBS curve ! Shapes procedure :: set_circle !!> Set a circle @@ -219,7 +223,7 @@ pure subroutine compute_dTgc_bspline_1d_vector(f_Xt, f_knot, f_degree, f_nc, f_n real(rk), allocatable, intent(out) :: f_Tgc(:,:) end subroutine - pure subroutine compute_dTgc_nurbs_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc) + pure subroutine compute_dTgc_nurbs_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc, f_elem) import :: rk real(rk), intent(in) :: f_Xt real(rk), intent(in), contiguous :: f_knot(:) @@ -228,9 +232,10 @@ pure subroutine compute_dTgc_nurbs_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_Wc, real(rk), intent(in), contiguous :: f_Wc(:) real(rk), allocatable, intent(out) :: f_dTgc(:) real(rk), allocatable, intent(out) :: f_Tgc(:) + integer, intent(in), optional :: f_elem(:) end subroutine - pure subroutine compute_dTgc_bspline_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_dTgc, f_Tgc) + pure subroutine compute_dTgc_bspline_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_dTgc, f_Tgc, f_elem) import :: rk real(rk), intent(in) :: f_Xt real(rk), intent(in), contiguous :: f_knot(:) @@ -238,6 +243,7 @@ pure subroutine compute_dTgc_bspline_1d_scalar(f_Xt, f_knot, f_degree, f_nc, f_d integer, intent(in) :: f_nc real(rk), allocatable, intent(out) :: f_dTgc(:) real(rk), allocatable, intent(out) :: f_Tgc(:) + integer, intent(in), optional :: f_elem(:) end subroutine end interface @@ -324,6 +330,36 @@ pure subroutine set1(this, knot, Xc, Wc) !=============================================================================== + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + !> Set knot vector, control points and weights for the NURBS curve object. + pure subroutine set1a(this, knot, Xc, Wc) + class(nurbs_curve), intent(inout) :: this + real(rk), intent(in), contiguous :: knot(:) + real(rk), intent(in), contiguous :: Xc(:) + real(rk), intent(in), contiguous, optional :: Wc(:) + + if (allocated(this%knot)) deallocate(this%knot) + if (allocated(this%Xc)) deallocate(this%Xc) + + this%knot = knot + call this%cmp_degree() + allocate(this%Xc(size(Xc), 3), source = 0.0_rk) + this%Xc(:,1) = Xc + this%nc = size(this%Xc, 1) + if (present(Wc)) then + if (size(Wc) /= this%nc) then + error stop 'Number of weights does not match the number of control points.' + else + if (allocated(this%Wc)) deallocate(this%Wc) + this%Wc = Wc + end if + end if + end subroutine + !=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -388,6 +424,44 @@ pure subroutine set3(this, Xc, Wc) !=============================================================================== + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine set4(this, degree, nc, Xc, Wc) + class(nurbs_curve), intent(inout) :: this + integer, intent(in) :: degree + integer, intent(in) :: nc + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), intent(in), contiguous, optional :: Wc(:) + integer :: m, i + + if (allocated(this%Xc)) deallocate(this%Xc) + + this%Xc = Xc + this%nc = nc + this%degree = degree + + ! Size of knot vectors + m = nc + degree + 1 + + if (allocated(this%knot)) deallocate(this%knot) + allocate(this%knot(m)) + this%knot(1:degree+1) = 0.0_rk + this%knot(degree+2:m-degree-1) = [(real(i, rk)/(m-2*degree-1), i=1, m-2*degree-2)] + this%knot(m-degree:m) = 1.0_rk + + if (present(Wc)) then + if (size(Wc) /= nc) then + error stop 'Number of weights does not match the number of control points.' + else + if (allocated(this%Wc)) deallocate(this%Wc) + this%Wc = Wc + end if + end if + end subroutine + !=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -1121,16 +1195,17 @@ pure subroutine derivative_vector(this, res, Xt, dTgc, Tgc) !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine derivative_scalar(this, Xt, dTgc, Tgc) + pure subroutine derivative_scalar(this, Xt, dTgc, Tgc, elem) class(nurbs_curve), intent(inout) :: this real(rk), intent(in) :: Xt + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:) real(rk), allocatable, intent(out), optional :: Tgc(:) if (this%is_rational()) then ! NURBS - call compute_dTgc(Xt, this%knot, this%degree, this%nc, this%Wc, dTgc, Tgc) + call compute_dTgc(Xt, this%knot, this%degree, this%nc, this%Wc, dTgc, Tgc, elem) else ! B-Spline - call compute_dTgc(Xt, this%knot, this%degree, this%nc, dTgc, Tgc) + call compute_dTgc(Xt, this%knot, this%degree, this%nc, dTgc, Tgc, elem) end if end subroutine !=============================================================================== @@ -1837,6 +1912,78 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest end subroutine !=============================================================================== + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine ansatz(this, ie, ig, Tgc, dTgc_dXg, dL) + class(nurbs_curve), intent(inout) :: this + + integer, intent(in) :: ie, ig + real(rk), intent(out) :: dL + real(rk), allocatable, intent(out) :: Tgc(:), dTgc_dXg(:,:) + real(rk), allocatable :: Xth(:), Xth_e(:), Xc_eT(:,:), Xksi(:), Wksi(:) + integer, allocatable :: elem_th(:,:), elem_c(:,:), elem_ce(:) + type(nurbs_curve) :: th, th_e + real(rk), allocatable :: dTtth_dXksi(:), Ttth(:), dTgc_dXt(:), dXg_dXt(:) + real(rk) :: Xt, dXt_dXksi + real(rk), allocatable :: dXg_dXksi(:) !> Jacobian matrix + real(rk) :: det_dXg_dXksi !> Determinant of the Jacobian matrix + + call gauss_leg([0.0_rk, 1.0_rk], this%degree, Xksi, Wksi) + + Xth = unique(this%knot) + + call th%set([0.0_rk,Xth,1.0_rk], Xth) + elem_th = th%cmp_elem() + + elem_c = this%cmp_elem() + + Xth_e = Xth(elem_th(ie,:)) + call th_e%set([0.0_rk,0.0_rk,1.0_rk,1.0_rk], Xth_e) + + elem_ce = elem_c(ie,:) + Xc_eT = transpose(this%Xc(elem_ce,:)) + + call th_e%derivative(Xksi(ig), dTtth_dXksi, Ttth) + Xt = dot_product(Xth_e, Ttth) + dXt_dXksi = dot_product(Xth_e, dTtth_dXksi) + + call this%derivative(Xt, dTgc_dXt, Tgc, elem_ce) + dXg_dXt = matmul(Xc_eT, dTgc_dXt) + + dTgc_dXg = dyad(dTgc_dXt, dXg_dXt)/norm2(dXg_dXt) + + dXg_dXksi = dXg_dXt*dXt_dXksi + det_dXg_dXksi = norm2(dXg_dXksi) + + dL = det_dXg_dXksi*Wksi(ig) + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine cmp_length(this, length) + class(nurbs_curve), intent(inout) :: this + real(rk), intent(out) :: length + real(rk), allocatable :: Tgc(:), dTgc_dXg(:,:) + integer :: ie, ig + real(rk) :: dL, dL_ig + + length = 0.0_rk + do ie = 1, size(this%cmp_elem(),1) + dL = 0.0_rk + do ig = 1, size(this%cmp_elem(),2) + call this%ansatz(ie, ig, Tgc, dTgc_dXg, dL_ig) + dL = dL + dL_ig + end do + length = length + dL + end do + end subroutine + !=============================================================================== + end module forcad_nurbs_curve !=============================================================================== @@ -1972,7 +2119,7 @@ impure subroutine compute_dTgc_nurbs_1d_vector(Xt, knot, degree, nc, ng, Wc, dTg !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure subroutine compute_dTgc_nurbs_1d_scalar(Xt, knot, degree, nc, Wc, dTgc, Tgc) +impure subroutine compute_dTgc_nurbs_1d_scalar(Xt, knot, degree, nc, Wc, dTgc, Tgc, elem) use forcad_utils, only: rk, basis_bspline_der implicit none @@ -1981,14 +2128,22 @@ impure subroutine compute_dTgc_nurbs_1d_scalar(Xt, knot, degree, nc, Wc, dTgc, T integer, intent(in) :: degree integer, intent(in) :: nc real(rk), intent(in), contiguous :: Wc(:) + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:) real(rk), allocatable, intent(out) :: Tgc(:) real(rk), allocatable :: dBi(:), Bi(:) - allocate(dTgc(nc), Tgc(nc)) call basis_bspline_der(Xt, knot, nc, degree, dBi, Bi) - Tgc = Bi*(Wc/(dot_product(Bi,Wc))) - dTgc = ( dBi*Wc - Tgc*dot_product(dBi,Wc) ) / dot_product(Bi,Wc) + + if (.not. present(elem)) then + allocate(dTgc(nc), Tgc(nc)) + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + dTgc = ( dBi*Wc - Tgc*dot_product(dBi,Wc) ) / dot_product(Bi,Wc) + else + allocate(dTgc(size(elem)), Tgc(size(elem))) + Tgc = Bi(elem)*(Wc(elem)/(dot_product(Bi(elem),Wc(elem)))) + dTgc = ( dBi(elem)*Wc(elem) - Tgc*dot_product(dBi(elem),Wc(elem)) ) / dot_product(Bi(elem),Wc(elem)) + end if end subroutine !=============================================================================== @@ -2023,7 +2178,7 @@ impure subroutine compute_dTgc_bspline_1d_vector(Xt, knot, degree, nc, ng, dTgc, !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure subroutine compute_dTgc_bspline_1d_scalar(Xt, knot, degree, nc, dTgc, Tgc) +impure subroutine compute_dTgc_bspline_1d_scalar(Xt, knot, degree, nc, dTgc, Tgc, elem) use forcad_utils, only: rk, basis_bspline_der implicit none @@ -2031,11 +2186,20 @@ impure subroutine compute_dTgc_bspline_1d_scalar(Xt, knot, degree, nc, dTgc, Tgc real(rk), intent(in), contiguous :: knot(:) integer, intent(in) :: degree integer, intent(in) :: nc + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:) real(rk), allocatable, intent(out) :: Tgc(:) - - allocate(dTgc(nc), Tgc(nc)) - call basis_bspline_der(Xt, knot, nc, degree, dTgc, Tgc) + real(rk), allocatable :: dB(:), B(:) + + if (.not. present(elem)) then + allocate(dTgc(nc), Tgc(nc)) + call basis_bspline_der(Xt, knot, nc, degree, dTgc, Tgc) + else + allocate(dB(size(elem)), B(size(elem))) + call basis_bspline_der(Xt, knot, nc, degree, dB, B) + Tgc = B(elem) + dTgc = dB(elem) + end if end subroutine !=============================================================================== diff --git a/src/forcad_nurbs_surface.f90 b/src/forcad_nurbs_surface.f90 index 9ab8a6328..4865dd9a4 100644 --- a/src/forcad_nurbs_surface.f90 +++ b/src/forcad_nurbs_surface.f90 @@ -5,7 +5,7 @@ module forcad_nurbs_surface use forcad_utils, only: rk, basis_bspline, elemConn_C0, kron, ndgrid, compute_multiplicity, compute_knot_vector, & basis_bspline_der, insert_knot_A_5_1, findspan, elevate_degree_A_5_9, remove_knots_A_5_8, tetragon_Xc, & - elemConn_Cn, unique, rotation + elemConn_Cn, unique, rotation, det, inv, gauss_leg implicit none @@ -34,7 +34,8 @@ module forcad_nurbs_surface procedure :: set1 !!> Set knot vectors, control points and weights for the NURBS surface object procedure :: set2 !!> Set NURBS surface using nodes of parameter space, degree, continuity, control points and weights procedure :: set3 !!> Set Bezier or Rational Bezier surface using control points and weights - generic :: set => set1, set2, set3 !!> Set NURBS surface + procedure :: set4 !!> Set NURBS surface using degree, number of control points, control points and weights + generic :: set => set1, set2, set3, set4 !!> Set NURBS surface procedure :: create !!> Generate geometry points procedure :: cmp_Xg !!> Compute geometry points procedure, private :: get_Xc_all !!> Get all control points @@ -97,6 +98,8 @@ module forcad_nurbs_surface procedure :: show !!> Show the NURBS object using PyVista procedure :: nearest_point !!> Find the nearest point on the NURBS surface (Approximation) procedure :: nearest_point2 !!> Find the nearest point on the NURBS surface (Minimization - Newton's method) + procedure :: ansatz !!> Compute the shape functions, derivative of shape functions and dA + procedure :: cmp_area !!> Compute the area of the NURBS surface ! Shapes procedure :: set_tetragon !!> Set a tetragon @@ -178,7 +181,7 @@ pure subroutine compute_dTgc_bspline_2d_vector(f_Xt, f_knot1, f_knot2, f_degree, real(rk), allocatable, intent(out) :: f_Tgc(:,:) end subroutine - pure subroutine compute_dTgc_nurbs_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc) + pure subroutine compute_dTgc_nurbs_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc, elem) import :: rk real(rk), intent(in), contiguous :: f_Xt(:) real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) @@ -187,9 +190,10 @@ pure subroutine compute_dTgc_nurbs_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, f real(rk), intent(in), contiguous :: f_Wc(:) real(rk), allocatable, intent(out) :: f_dTgc(:,:) real(rk), allocatable, intent(out) :: f_Tgc(:) + integer, intent(in) :: elem(:) end subroutine - pure subroutine compute_dTgc_bspline_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, nc, f_dTgc, f_Tgc) + pure subroutine compute_dTgc_bspline_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, nc, f_dTgc, f_Tgc, elem) import :: rk real(rk), intent(in), contiguous :: f_Xt(:) real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:) @@ -197,6 +201,7 @@ pure subroutine compute_dTgc_bspline_2d_scalar(f_Xt, f_knot1, f_knot2, f_degree, integer, intent(in) :: nc(2) real(rk), allocatable, intent(out) :: f_dTgc(:,:) real(rk), allocatable, intent(out) :: f_Tgc(:) + integer, intent(in) :: elem(:) end subroutine end interface @@ -397,6 +402,50 @@ pure subroutine set3(this, nc, Xc, Wc) !=============================================================================== + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine set4(this, degree, nc, Xc, Wc) + class(nurbs_surface), intent(inout) :: this + integer, intent(in), contiguous :: degree(:) + integer, intent(in), contiguous :: nc(:) + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), intent(in), contiguous, optional :: Wc(:) + integer :: m(3), i + + if (allocated(this%Xc)) deallocate(this%Xc) + + this%Xc = Xc + this%nc = nc + this%degree = degree + + ! Size of knot vectors + m = nc + degree + 1 + + if (allocated(this%knot1)) deallocate(this%knot1) + allocate(this%knot1(m(1))) + this%knot1(1:degree(1)+1) = 0.0_rk + this%knot1(degree(1)+2:m(1)-degree(1)-1) = [(real(i, rk)/(m(1)-2*degree(1)-1), i=1, m(1)-2*degree(1)-2)] + this%knot1(m(1)-degree(1):m(1)) = 1.0_rk + + if (allocated(this%knot2)) deallocate(this%knot2) + allocate(this%knot2(m(2))) + this%knot2(1:degree(2)+1) = 0.0_rk + this%knot2(degree(2)+2:m(2)-degree(2)-1) = [(real(i, rk)/(m(2)-2*degree(2)-1), i=1, m(2)-2*degree(2)-2)] + this%knot2(m(2)-degree(2):m(2)) = 1.0_rk + + if (present(Wc)) then + if (size(Wc) /= nc(1)*nc(2)) then + error stop 'Number of weights does not match the number of control points.' + else + if (allocated(this%Wc)) deallocate(this%Wc) + this%Wc = Wc + end if + end if + end subroutine + !=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -1228,16 +1277,17 @@ pure subroutine derivative_vector(this, res1, res2, Xt1, Xt2, dTgc, Tgc) !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine derivative_scalar(this, Xt, dTgc, Tgc) + pure subroutine derivative_scalar(this, Xt, dTgc, Tgc, elem) class(nurbs_surface), intent(inout) :: this real(rk), intent(in), contiguous :: Xt(:) + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:,:) real(rk), allocatable, intent(out), optional :: Tgc(:) if (this%is_rational()) then ! NURBS - call compute_dTgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%Wc, dTgc, Tgc) + call compute_dTgc(Xt, this%knot1, this%knot2, this%degree, this%nc, this%Wc, dTgc, Tgc, elem) else ! B-Spline - call compute_dTgc(Xt, this%knot1, this%knot2, this%degree, this%nc, dTgc, Tgc) + call compute_dTgc(Xt, this%knot1, this%knot2, this%degree, this%nc, dTgc, Tgc, elem) end if end subroutine !=============================================================================== @@ -2395,7 +2445,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest integer, intent(in) :: maxit real(rk), intent(out) :: nearest_Xt(2) real(rk), allocatable, intent(out), optional :: nearest_Xg(:) - real(rk):: obj, grad(2), hess(2,2), dk(2), alphak, tau, beta, det_inv, Ainv(2,2), lower_bounds(2), upper_bounds(2) + real(rk):: obj, grad(2), hess(2,2), dk(2), alphak, tau, beta, lower_bounds(2), upper_bounds(2) real(rk), allocatable :: Xg(:), xk(:), Tgc(:), dTgc(:,:), d2Tgc(:,:) integer :: k, l logical :: convergenz @@ -2465,15 +2515,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt) else - ! Inverse of Hessian - det_inv = 1.0_rk/(hess(1,1)*hess(2,2) - hess(1,2)*hess(2,1)) - Ainv(1,1) = hess(2,2) - Ainv(2,1) = -hess(2,1) - Ainv(1,2) = -hess(1,2) - Ainv(2,2) = hess(1,1) - Ainv = Ainv * det_inv - - dk = - matmul(Ainv, grad) + dk = - matmul(inv(hess), grad) ! Backtracking-Armijo Line Search alphak = 1.0_rk @@ -2495,6 +2537,79 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest end subroutine !=============================================================================== + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine ansatz(this, ie, ig, Tgc, dTgc_dXg, dA) + class(nurbs_surface), intent(inout) :: this + integer, intent(in) :: ie, ig + real(rk), intent(out) :: dA + real(rk), allocatable, intent(out) :: Tgc(:), dTgc_dXg(:,:) + real(rk), allocatable :: Xth(:,:), Xth_e(:,:), Xth_eT(:,:), Xc_eT(:,:), Xth1(:), Xth2(:), Xksi(:,:), Wksi(:) + integer, allocatable :: elem_th(:,:), elem_c(:,:), elem_ce(:) + type(nurbs_surface) :: th, th_e + real(rk), allocatable :: dTtth_dXksi(:,:), Ttth(:), dTgc_dXt(:,:), Xt(:), dXt_dXksi(:,:), dXg_dXt(:,:) + real(rk), allocatable :: dXg_dXksi(:,:) !> Jacobian matrix + real(rk) :: det_dXg_dXksi !> Determinant of the Jacobian matrix + + call gauss_leg([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], this%degree, Xksi, Wksi) + + Xth1 = unique(this%knot1) + Xth2 = unique(this%knot2) + call ndgrid(Xth1, Xth2, Xth) + + call th%set([0.0_rk,Xth1,1.0_rk], [0.0_rk,Xth2,1.0_rk], Xth) + elem_th = th%cmp_elem() + + elem_c = this%cmp_elem() + + Xth_e = Xth(elem_th(ie,:),:) + call th_e%set([0.0_rk,0.0_rk,1.0_rk,1.0_rk], [0.0_rk,0.0_rk,1.0_rk,1.0_rk], Xth_e) + + Xth_eT = transpose(Xth_e) + elem_ce = elem_c(ie,:) + Xc_eT = transpose(this%Xc(elem_ce,:)) + + call th_e%derivative(Xksi(ig,:), dTtth_dXksi, Ttth) + Xt = matmul(Xth_eT, Ttth) + dXt_dXksi = matmul(Xth_eT, dTtth_dXksi) + + call this%derivative(Xt, dTgc_dXt, Tgc, elem_ce) + dXg_dXt = matmul(Xc_eT, dTgc_dXt) + + dTgc_dXg = matmul(dTgc_dXt, inv(dXg_dXt)) + + dXg_dXksi = matmul(dXg_dXt, dXt_dXksi) + det_dXg_dXksi = det(dXg_dXksi) + + dA = det_dXg_dXksi*Wksi(ig) + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine cmp_area(this, area) + class(nurbs_surface), intent(inout) :: this + real(rk), intent(out) :: area + real(rk), allocatable :: Tgc(:), dTgc_dXg(:,:) + integer :: ie, ig + real(rk) :: dA, dA_ig + + area = 0.0_rk + do ie = 1, size(this%cmp_elem(),1) + dA = 0.0_rk + do ig = 1, size(this%cmp_elem(),2) + call this%ansatz(ie, ig, Tgc, dTgc_dXg, dA_ig) + dA = dA + dA_ig + end do + area = area + dA + end do + end subroutine + !=============================================================================== + end module forcad_nurbs_surface !=============================================================================== @@ -2653,7 +2768,7 @@ impure subroutine compute_dTgc_nurbs_2d_vector(Xt, knot1, knot2, degree, nc, ng, !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure subroutine compute_dTgc_nurbs_2d_scalar(Xt, knot1, knot2, degree, nc, Wc, dTgc, Tgc) +impure subroutine compute_dTgc_nurbs_2d_scalar(Xt, knot1, knot2, degree, nc, Wc, dTgc, Tgc, elem) use forcad_utils, only: rk, basis_bspline_der, kron implicit none @@ -2662,25 +2777,44 @@ impure subroutine compute_dTgc_nurbs_2d_scalar(Xt, knot1, knot2, degree, nc, Wc, integer, intent(in) :: degree(2) integer, intent(in) :: nc(2) real(rk), intent(in), contiguous :: Wc(:) + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:,:) real(rk), allocatable, intent(out) :: Tgc(:) real(rk), allocatable :: dB1(:), dB2(:), dBi(:,:) real(rk), allocatable :: B1(:), B2(:), Bi(:) - allocate(dTgc(nc(1)*nc(2), 2), Tgc(nc(1)*nc(2))) - allocate(dBi(nc(1)*nc(2), 2), Bi(nc(1)*nc(2))) - call basis_bspline_der(Xt(1), knot1, nc(1), degree(1), dB1, B1) call basis_bspline_der(Xt(2), knot2, nc(2), degree(2), dB2, B2) - Bi = kron(B2, B1) - Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + if (.not. present(elem)) then + allocate(dTgc(nc(1)*nc(2), 2), Tgc(nc(1)*nc(2))) + allocate(dBi(nc(1)*nc(2), 2), Bi(nc(1)*nc(2))) - dBi(:,1) = kron(B2, dB1) - dBi(:,2) = kron(dB2, B1) + Bi = kron(B2, B1) + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) - dTgc(:,1) = ( dBi(:,1)*Wc - Tgc*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) - dTgc(:,2) = ( dBi(:,2)*Wc - Tgc*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + dBi(:,1) = kron(B2, dB1) + dBi(:,2) = kron(dB2, B1) + + dTgc(:,1) = ( dBi(:,1)*Wc - Tgc*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(:,2) = ( dBi(:,2)*Wc - Tgc*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + else + allocate(dTgc(size(elem), 2), Tgc(size(elem))) + allocate(dBi(size(elem), 2), Bi(size(elem))) + + associate(Biall => kron(B2, B1)) + Bi = Biall(elem) + Tgc = Bi*(Wc(elem)/(dot_product(Bi,Wc(elem)))) + end associate + + associate(dB1all => kron(B2, dB1), dB2all => kron(dB2, B1)) + dBi(:,1) = dB1all(elem) + dBi(:,2) = dB2all(elem) + end associate + + dTgc(:,1) = ( dBi(:,1)*Wc(elem) - Tgc*dot_product(dBi(:,1),Wc(elem)) ) / dot_product(Bi,Wc(elem)) + dTgc(:,2) = ( dBi(:,2)*Wc(elem) - Tgc*dot_product(dBi(:,2),Wc(elem)) ) / dot_product(Bi,Wc(elem)) + end if end subroutine !=============================================================================== @@ -2721,7 +2855,7 @@ impure subroutine compute_dTgc_bspline_2d_vector(Xt, knot1, knot2, degree, nc, n !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure subroutine compute_dTgc_bspline_2d_scalar(Xt, knot1, knot2, degree, nc, dTgc, Tgc) +impure subroutine compute_dTgc_bspline_2d_scalar(Xt, knot1, knot2, degree, nc, dTgc, Tgc, elem) use forcad_utils, only: rk, basis_bspline_der, kron implicit none @@ -2729,18 +2863,33 @@ impure subroutine compute_dTgc_bspline_2d_scalar(Xt, knot1, knot2, degree, nc, d real(rk), intent(in), contiguous :: knot1(:), knot2(:) integer, intent(in) :: degree(2) integer, intent(in) :: nc(2) + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:,:) real(rk), allocatable, intent(out) :: Tgc(:) real(rk), allocatable :: dTgc1(:), dTgc2(:) real(rk), allocatable :: Tgc1(:), Tgc2(:) - allocate(dTgc(nc(1)*nc(2), 2)) call basis_bspline_der(Xt(1), knot1, nc(1), degree(1), dTgc1, Tgc1) call basis_bspline_der(Xt(2), knot2, nc(2), degree(2), dTgc2, Tgc2) - Tgc = kron(Tgc2, Tgc1) - dTgc(:,1) = kron(Tgc2, dTgc1) - dTgc(:,2) = kron(dTgc2, Tgc1) + if (.not. present(elem)) then + allocate(dTgc(nc(1)*nc(2), 2)) + Tgc = kron(Tgc2, Tgc1) + + dTgc(:,1) = kron(Tgc2, dTgc1) + dTgc(:,2) = kron(dTgc2, Tgc1) + else + allocate(dTgc(size(elem), 2)) + + associate(B => kron(Tgc2, Tgc1)) + Tgc = B(elem) + end associate + + associate(dB1 => kron(Tgc2, dTgc1), dB2 => kron(dTgc2, Tgc1)) + dTgc(:,1) = dB1(elem) + dTgc(:,2) = dB2(elem) + end associate + end if end subroutine !=============================================================================== diff --git a/src/forcad_nurbs_volume.f90 b/src/forcad_nurbs_volume.f90 index b156f08b3..97eeae690 100644 --- a/src/forcad_nurbs_volume.f90 +++ b/src/forcad_nurbs_volume.f90 @@ -5,7 +5,7 @@ module forcad_nurbs_volume use forcad_utils, only: rk, basis_bspline, elemConn_C0, kron, ndgrid, compute_multiplicity, compute_knot_vector, & basis_bspline_der, insert_knot_A_5_1, findspan, elevate_degree_A_5_9, hexahedron_Xc, remove_knots_A_5_8, & - elemConn_Cn, unique, rotation + elemConn_Cn, unique, rotation, det, inv, gauss_leg implicit none @@ -36,7 +36,8 @@ module forcad_nurbs_volume procedure :: set1 !!> Set knot vectors, control points and weights for the NURBS volume object procedure :: set2 !!> Set NURBS volume using nodes of parameter space, degree, continuity, control points and weights procedure :: set3 !!> Set Bezier or Rational Bezier volume using control points and weights - generic :: set => set1, set2, set3 !!> Set NURBS volume + procedure :: set4 !!> Set NURBS volume using degree, number of control points, control points and weights + generic :: set => set1, set2, set3, set4 !!> Set NURBS volume procedure :: create !!> Generate geometry points procedure :: cmp_Xg !!> Compute geometry points procedure, private :: get_Xc_all !!> Get all control points @@ -100,6 +101,9 @@ module forcad_nurbs_volume procedure :: show !!> Show the NURBS object using PyVista procedure :: nearest_point !!> Find the nearest point on the NURBS volume (Approximation) procedure :: nearest_point2 !!> Find the nearest point on the NURBS volume (Minimization - Newton's method) + procedure :: ansatz !!> Compute the shape functions, derivative of shape functions and dV + procedure :: cmp_volume !!> Compute the volume of the NURBS volume + ! Faces procedure :: cmp_elemFace_Xc_vis !!> Compute faces of the control points @@ -185,7 +189,7 @@ pure subroutine compute_dTgc_bspline_3d_vector(f_Xt, f_knot1, f_knot2, f_knot3, real(rk), allocatable, intent(out) :: f_Tgc(:,:) end subroutine - pure subroutine compute_dTgc_nurbs_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc) + pure subroutine compute_dTgc_nurbs_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_Wc, f_dTgc, f_Tgc, f_elem) import :: rk real(rk), intent(in), contiguous :: f_Xt(:) real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) @@ -194,9 +198,10 @@ pure subroutine compute_dTgc_nurbs_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_ real(rk), intent(in), contiguous :: f_Wc(:) real(rk), allocatable, intent(out) :: f_dTgc(:,:) real(rk), allocatable, intent(out) :: f_Tgc(:) + integer, intent(in), optional :: f_elem(:) end subroutine - pure subroutine compute_dTgc_bspline_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_dTgc, f_Tgc) + pure subroutine compute_dTgc_bspline_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f_nc, f_dTgc, f_Tgc, f_elem) import :: rk real(rk), intent(in), contiguous :: f_Xt(:) real(rk), intent(in), contiguous :: f_knot1(:), f_knot2(:), f_knot3(:) @@ -204,6 +209,7 @@ pure subroutine compute_dTgc_bspline_3d_scalar(f_Xt, f_knot1, f_knot2, f_knot3, integer, intent(in) :: f_nc(3) real(rk), allocatable, intent(out) :: f_dTgc(:,:) real(rk), allocatable, intent(out) :: f_Tgc(:) + integer, intent(in), optional :: f_elem(:) end subroutine end interface @@ -414,6 +420,56 @@ pure subroutine set3(this, nc, Xc, Wc) !=============================================================================== + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine set4(this, degree, nc, Xc, Wc) + class(nurbs_volume), intent(inout) :: this + integer, intent(in), contiguous :: degree(:) + integer, intent(in), contiguous :: nc(:) + real(rk), intent(in), contiguous :: Xc(:,:) + real(rk), intent(in), contiguous, optional :: Wc(:) + integer :: m(3), i + + if (allocated(this%Xc)) deallocate(this%Xc) + + this%Xc = Xc + this%nc = nc + this%degree = degree + + ! Size of knot vectors + m = nc + degree + 1 + + if (allocated(this%knot1)) deallocate(this%knot1) + allocate(this%knot1(m(1))) + this%knot1(1:degree(1)+1) = 0.0_rk + this%knot1(degree(1)+2:m(1)-degree(1)-1) = [(real(i, rk)/(m(1)-2*degree(1)-1), i=1, m(1)-2*degree(1)-2)] + this%knot1(m(1)-degree(1):m(1)) = 1.0_rk + + if (allocated(this%knot2)) deallocate(this%knot2) + allocate(this%knot2(m(2))) + this%knot2(1:degree(2)+1) = 0.0_rk + this%knot2(degree(2)+2:m(2)-degree(2)-1) = [(real(i, rk)/(m(2)-2*degree(2)-1), i=1, m(2)-2*degree(2)-2)] + this%knot2(m(2)-degree(2):m(2)) = 1.0_rk + + if (allocated(this%knot3)) deallocate(this%knot3) + allocate(this%knot3(m(3))) + this%knot3(1:degree(3)+1) = 0.0_rk + this%knot3(degree(3)+2:m(3)-degree(3)-1) = [(real(i, rk)/(m(3)-2*degree(3)-1), i=1, m(3)-2*degree(3)-2)] + this%knot3(m(3)-degree(3):m(3)) = 1.0_rk + + if (present(Wc)) then + if (size(Wc) /= nc(1)*nc(2)*nc(3)) then + error stop 'Number of weights does not match the number of control points.' + else + if (allocated(this%Wc)) deallocate(this%Wc) + this%Wc = Wc + end if + end if + end subroutine + !=============================================================================== + + !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause @@ -1354,16 +1410,17 @@ pure subroutine derivative_vector(this, res1, res2, res3, Xt1, Xt2, Xt3, dTgc, T !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause - pure subroutine derivative_scalar(this, Xt, dTgc, Tgc) + pure subroutine derivative_scalar(this, Xt, dTgc, Tgc, elem) class(nurbs_volume), intent(inout) :: this real(rk), intent(in), contiguous :: Xt(:) + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:,:) real(rk), allocatable, intent(out), optional :: Tgc(:) if (this%is_rational()) then ! NURBS - call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%Wc, dTgc, Tgc) + call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%Wc, dTgc, Tgc, elem) else - call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, dTgc, Tgc) + call compute_dTgc(Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, dTgc, Tgc, elem) end if end subroutine !=============================================================================== @@ -2909,7 +2966,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest integer, intent(in) :: maxit real(rk), intent(out) :: nearest_Xt(3) real(rk), allocatable, intent(out), optional :: nearest_Xg(:) - real(rk):: obj, grad(3), hess(3,3), dk(3), alphak, tau, beta, det_inv, Ainv(3,3), lower_bounds(3), upper_bounds(3) + real(rk):: obj, grad(3), hess(3,3), dk(3), alphak, tau, beta, lower_bounds(3), upper_bounds(3) real(rk), allocatable :: Xg(:), xk(:), Tgc(:), dTgc(:,:), d2Tgc(:,:) integer :: k, l logical :: convergenz @@ -3001,24 +3058,7 @@ impure subroutine nearest_point2(this, point_Xg, tol, maxit, nearest_Xt, nearest if (present(nearest_Xg)) nearest_Xg = this%cmp_Xg(nearest_Xt) else - ! Inverse of Hessian - det_inv = 1.0_rk/(& - + hess(1,1)*hess(2,2)*hess(3,3) - hess(1,1)*hess(2,3)*hess(3,2)& - - hess(1,2)*hess(2,1)*hess(3,3) + hess(1,2)*hess(2,3)*hess(3,1)& - + hess(1,3)*hess(2,1)*hess(3,2) - hess(1,3)*hess(2,2)*hess(3,1)) - Ainv(1,1) = +(hess(2,2)*hess(3,3) - hess(2,3)*hess(3,2)) - Ainv(2,1) = -(hess(2,1)*hess(3,3) - hess(2,3)*hess(3,1)) - Ainv(3,1) = +(hess(2,1)*hess(3,2) - hess(2,2)*hess(3,1)) - Ainv(1,2) = -(hess(1,2)*hess(3,3) - hess(1,3)*hess(3,2)) - Ainv(2,2) = +(hess(1,1)*hess(3,3) - hess(1,3)*hess(3,1)) - Ainv(3,2) = -(hess(1,1)*hess(3,2) - hess(1,2)*hess(3,1)) - Ainv(1,3) = +(hess(1,2)*hess(2,3) - hess(1,3)*hess(2,2)) - Ainv(2,3) = -(hess(1,1)*hess(2,3) - hess(1,3)*hess(2,1)) - Ainv(3,3) = +(hess(1,1)*hess(2,2) - hess(1,2)*hess(2,1)) - - Ainv = det_inv*Ainv - - dk = - matmul(Ainv, grad) + dk = - matmul(inv(hess), grad) ! Backtracking-Armijo Line Search alphak = 1.0_rk @@ -3227,6 +3267,79 @@ pure function cmp_degreeFace(this, face) result(degree) end function !=============================================================================== + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine ansatz(this, ie, ig, Tgc, dTgc_dXg, dV) + class(nurbs_volume), intent(inout) :: this + integer, intent(in) :: ie, ig + real(rk), intent(out) :: dV + real(rk), allocatable, intent(out) :: Tgc(:), dTgc_dXg(:,:) + real(rk), allocatable :: Xth(:,:), Xth_e(:,:), Xth_eT(:,:), Xc_eT(:,:), Xth1(:), Xth2(:), Xth3(:), Xksi(:,:), Wksi(:) + integer, allocatable :: elem_th(:,:), elem_c(:,:), elem_ce(:) + type(nurbs_volume) :: th, th_e + real(rk), allocatable :: dTtth_dXksi(:,:), Ttth(:), dTgc_dXt(:,:), Xt(:), dXt_dXksi(:,:), dXg_dXt(:,:) + real(rk), allocatable :: dXg_dXksi(:,:) !> Jacobian matrix + real(rk) :: det_dXg_dXksi !> Determinant of the Jacobian matrix + + call gauss_leg([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], this%degree, Xksi, Wksi) + + Xth1 = unique(this%knot1) + Xth2 = unique(this%knot2) + Xth3 = unique(this%knot3) + call ndgrid(Xth1, Xth2, Xth3, Xth) + call th%set([0.0_rk,Xth1,1.0_rk], [0.0_rk,Xth2,1.0_rk], [0.0_rk,Xth3,1.0_rk], Xth) + elem_th = th%cmp_elem() + + elem_c = this%cmp_elem() + + Xth_e = Xth(elem_th(ie,:),:) + call th_e%set([0.0_rk,0.0_rk,1.0_rk,1.0_rk], [0.0_rk,0.0_rk,1.0_rk,1.0_rk], [0.0_rk,0.0_rk,1.0_rk,1.0_rk], Xth_e) + + Xth_eT = transpose(Xth_e) + elem_ce = elem_c(ie,:) + Xc_eT = transpose(this%Xc(elem_ce,:)) + + call th_e%derivative(Xksi(ig,:), dTtth_dXksi, Ttth) + Xt = matmul(Xth_eT, Ttth) + dXt_dXksi = matmul(Xth_eT, dTtth_dXksi) + + call this%derivative(Xt, dTgc_dXt, Tgc, elem_ce) + dXg_dXt = matmul(Xc_eT, dTgc_dXt) + + dTgc_dXg = matmul(dTgc_dXt, inv(dXg_dXt)) + + dXg_dXksi = matmul(dXg_dXt, dXt_dXksi) + det_dXg_dXksi = det(dXg_dXksi) + + dV = det_dXg_dXksi*Wksi(ig) + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine cmp_volume(this, volume) + class(nurbs_volume), intent(inout) :: this + real(rk), intent(out) :: volume + real(rk), allocatable :: Tgc(:), dTgc_dXg(:,:) + integer :: ie, ig + real(rk) :: dV, dV_ig + + volume = 0.0_rk + do ie = 1, size(this%cmp_elem(),1) + dV = 0.0_rk + do ig = 1, size(this%cmp_elem(),2) + call this%ansatz(ie, ig, Tgc, dTgc_dXg, dV_ig) + dV = dV + dV_ig + end do + volume = volume + dV + end do + end subroutine + !=============================================================================== + end module forcad_nurbs_volume !=============================================================================== @@ -3390,7 +3503,7 @@ impure subroutine compute_dTgc_nurbs_3d_vector(Xt, knot1, knot2, knot3, degree, !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure subroutine compute_dTgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, Wc, dTgc, Tgc) +impure subroutine compute_dTgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, Wc, dTgc, Tgc, elem) use forcad_utils, only: rk, basis_bspline_der, kron implicit none @@ -3399,29 +3512,51 @@ impure subroutine compute_dTgc_nurbs_3d_scalar(Xt, knot1, knot2, knot3, degree, integer, intent(in) :: degree(3) integer, intent(in) :: nc(3) real(rk), intent(in), contiguous :: Wc(:) + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:,:) real(rk), allocatable, intent(out) :: Tgc(:) real(rk), allocatable :: dB1(:), dB2(:), dB3(:), dBi(:,:) real(rk), allocatable :: B1(:), B2(:), B3(:), Bi(:) - allocate(dTgc(nc(1)*nc(2)*nc(3), 3)) - allocate(Tgc(nc(1)*nc(2)*nc(3))) - allocate(dBi(nc(1)*nc(2)*nc(3), 3), Bi(nc(1)*nc(2)*nc(3))) - call basis_bspline_der(Xt(1), knot1, nc(1), degree(1), dB1, B1) call basis_bspline_der(Xt(2), knot2, nc(2), degree(2), dB2, B2) call basis_bspline_der(Xt(3), knot3, nc(3), degree(3), dB3, B3) - Bi = kron( B3, kron( B2, B1)) - Tgc = Bi*(Wc/(dot_product(Bi,Wc))) + if (.not. present(elem)) then + allocate(dTgc(nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(nc(1)*nc(2)*nc(3))) + allocate(dBi(nc(1)*nc(2)*nc(3), 3), Bi(nc(1)*nc(2)*nc(3))) - dBi(:,1) = kron(kron(B3,B2),dB1) - dBi(:,2) = kron(kron(B3,dB2),B1) - dBi(:,3) = kron(kron(dB3,B2),B1) + Bi = kron( B3, kron( B2, B1)) + Tgc = Bi*(Wc/(dot_product(Bi,Wc))) - dTgc(:,1) = ( dBi(:,1)*Wc - Tgc*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) - dTgc(:,2) = ( dBi(:,2)*Wc - Tgc*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) - dTgc(:,3) = ( dBi(:,3)*Wc - Tgc*dot_product(dBi(:,3),Wc) ) / dot_product(Bi,Wc) + dBi(:,1) = kron(kron(B3,B2),dB1) + dBi(:,2) = kron(kron(B3,dB2),B1) + dBi(:,3) = kron(kron(dB3,B2),B1) + + dTgc(:,1) = ( dBi(:,1)*Wc - Tgc*dot_product(dBi(:,1),Wc) ) / dot_product(Bi,Wc) + dTgc(:,2) = ( dBi(:,2)*Wc - Tgc*dot_product(dBi(:,2),Wc) ) / dot_product(Bi,Wc) + dTgc(:,3) = ( dBi(:,3)*Wc - Tgc*dot_product(dBi(:,3),Wc) ) / dot_product(Bi,Wc) + else + allocate(dTgc(size(elem), 3)) + allocate(Tgc(size(elem))) + allocate(dBi(size(elem), 3), Bi(size(elem))) + + associate(Biall => kron( B3, kron( B2, B1))) + Bi = Biall(elem) + Tgc = Bi*(Wc(elem)/(dot_product(Bi,Wc(elem)))) + end associate + + associate(dB1all => kron(kron(B3,B2),dB1), dB2all => kron(kron(B3,dB2),B1), dB3all => kron(kron(dB3,B2),B1)) + dBi(:,1) = dB1all(elem) + dBi(:,2) = dB2all(elem) + dBi(:,3) = dB3all(elem) + end associate + + dTgc(:,1) = ( dBi(:,1)*Wc(elem) - Tgc*dot_product(dBi(:,1),Wc(elem)) ) / dot_product(Bi,Wc(elem)) + dTgc(:,2) = ( dBi(:,2)*Wc(elem) - Tgc*dot_product(dBi(:,2),Wc(elem)) ) / dot_product(Bi,Wc(elem)) + dTgc(:,3) = ( dBi(:,3)*Wc(elem) - Tgc*dot_product(dBi(:,3),Wc(elem)) ) / dot_product(Bi,Wc(elem)) + end if end subroutine !=============================================================================== @@ -3466,7 +3601,7 @@ impure subroutine compute_dTgc_bspline_3d_vector(Xt, knot1, knot2, knot3, degree !=============================================================================== !> author: Seyed Ali Ghasemi !> license: BSD 3-Clause -impure subroutine compute_dTgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, dTgc, Tgc) +impure subroutine compute_dTgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree, nc, dTgc, Tgc, elem) use forcad_utils, only: rk, basis_bspline_der, kron implicit none @@ -3474,21 +3609,38 @@ impure subroutine compute_dTgc_bspline_3d_scalar(Xt, knot1, knot2, knot3, degree real(rk), intent(in), contiguous :: knot1(:), knot2(:), knot3(:) integer, intent(in) :: degree(3) integer, intent(in) :: nc(3) + integer, intent(in), optional :: elem(:) real(rk), allocatable, intent(out) :: dTgc(:,:) real(rk), allocatable, intent(out) :: Tgc(:) real(rk), allocatable :: dB1(:), dB2(:), dB3(:) real(rk), allocatable :: B1(:), B2(:), B3(:) - allocate(dTgc(nc(1)*nc(2)*nc(3), 3)) - allocate(Tgc(nc(1)*nc(2)*nc(3))) call basis_bspline_der(Xt(1), knot1, nc(1), degree(1), dB1, B1) call basis_bspline_der(Xt(2), knot2, nc(2), degree(2), dB2, B2) call basis_bspline_der(Xt(3), knot3, nc(3), degree(3), dB3, B3) - Tgc = kron(B3, kron(B2, B1)) - dTgc(:,1) = kron(kron(B3,B2),dB1) - dTgc(:,2) = kron(kron(B3,dB2),B1) - dTgc(:,3) = kron(kron(dB3,B2),B1) + if (.not. present(elem)) then + allocate(dTgc(nc(1)*nc(2)*nc(3), 3)) + allocate(Tgc(nc(1)*nc(2)*nc(3))) + + Tgc = kron(B3, kron(B2, B1)) + dTgc(:,1) = kron(kron(B3,B2),dB1) + dTgc(:,2) = kron(kron(B3,dB2),B1) + dTgc(:,3) = kron(kron(dB3,B2),B1) + else + allocate(dTgc(size(elem), 3)) + allocate(Tgc(size(elem))) + + associate(B => kron(B3, kron(B2, B1))) + Tgc = B(elem) + end associate + + associate(dB1 => kron(kron(B3,B2),dB1), dB2 => kron(kron(B3,dB2),B1), dB3 => kron(kron(dB3,B2),B1)) + dTgc(:,1) = dB1(elem) + dTgc(:,2) = dB2(elem) + dTgc(:,3) = dB3(elem) + end associate + end if end subroutine !=============================================================================== diff --git a/src/forcad_utils.f90 b/src/forcad_utils.f90 index 95cac9104..380d3227c 100644 --- a/src/forcad_utils.f90 +++ b/src/forcad_utils.f90 @@ -3,12 +3,14 @@ !> This module contains parameters, functions and subroutines that are used in the library. module forcad_utils + use stdlib_quadrature, only: gauss_legendre + implicit none private public :: rk, basis_bernstein, basis_bspline, elemConn_C0, kron, ndgrid, compute_multiplicity, compute_knot_vector, & basis_bspline_der, insert_knot_A_5_1, findspan, elevate_degree_A_5_9, hexahedron_Xc, tetragon_Xc, remove_knots_A_5_8, & - elemConn_Cn, unique, rotation, basis_bspline_2der + elemConn_Cn, unique, rotation, basis_bspline_2der, det, inv, dyad, gauss_leg integer, parameter :: rk = kind(1.0d0) @@ -53,6 +55,22 @@ module forcad_utils end interface !=============================================================================== + + !=============================================================================== + interface dyad + module procedure dyad_t1_t1 + end interface + !=============================================================================== + + + !=============================================================================== + interface gauss_leg + module procedure gauss_legendre_1D + module procedure gauss_legendre_2D + module procedure gauss_legendre_3D + end interface + !=============================================================================== + contains !=============================================================================== @@ -119,7 +137,7 @@ pure subroutine basis_bspline_der(Xt, knot, nc, degree, dB, B) end if if ( Xth_ip /= Xth_i ) then N(i,p) = (Xt - Xth_i)/(Xth_ip - Xth_i) * N(i,p-1) - dN_dXt(i,p) = ( N(i,p-1) + (Xt - Xth_i)*dN_dXt(i,p-1) ) / (Xth_ip - Xth_i) + dN_dXt(i,p) = ( N(i,p-1) + (Xt - Xth_i)*dN_dXt(i,p-1) ) / (Xth_ip - Xth_i) end if if ( Xth_ip1 /= Xth_i1 ) then N(i,p) = N(i,p) + (Xth_ip1 - Xt)/(Xth_ip1 - Xth_i1) * N(i+1,p-1) @@ -170,7 +188,7 @@ pure subroutine basis_bspline_2der(Xt, knot, nc, degree, d2B, dB, B) if ( Xth_ip /= Xth_i) then N(i,p) = (Xt - Xth_i)/(Xth_ip - Xth_i) * N(i,p-1) ! dN_dXt(i,p) = p*(N(i,p-1)/(Xth_ip - Xth_i)) - dN_dXt(i,p) = ( N(i,p-1) + (Xt - Xth_i)*dN_dXt(i,p-1) ) / (Xth_ip - Xth_i) + dN_dXt(i,p) = ( N(i,p-1) + (Xt - Xth_i)*dN_dXt(i,p-1) ) / (Xth_ip - Xth_i) d2N_dXt2(i,p) = (2*dN_dXt(i,p-1) + (Xt - Xth_i)*d2N_dXt2(i,p-1)) / (Xth_ip - Xth_i) end if if ( Xth_ip1 /= Xth_i1 ) then @@ -1080,4 +1098,150 @@ pure function rotation(alpha, beta, theta) result(R) end function !=============================================================================== + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function det(A) result(detA) + real(rk), intent(in) :: A(:,:) + real(rk) :: detA + + if (size(A,1) == size(A,2)) then + select case(size(A,1)) + case(2) + detA = A(1,1)*A(2,2) - A(1,2)*A(2,1) + case(3) + detA = & + + A(1,1)*( A(2,2)*A(3,3) - A(2,3)*A(3,2) )& + - A(1,2)*( A(2,1)*A(3,3) - A(2,3)*A(3,1) )& + + A(1,3)*( A(2,1)*A(3,2) - A(2,2)*A(3,1) ) + end select + elseif (size(A,1) == 3 .and. size(A,2) == 2) then + detA = & + + A(1,1) * ( A(2,2) * 1.0_rk - A(3,2) * 1.0_rk )& + - A(1,2) * ( A(2,1) * 1.0_rk - A(3,1) * 1.0_rk )& + + 1.0_rk * ( A(2,1) * A(3,2) - A(3,1) * A(2,2) ) + end if + end function + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + recursive pure function inv(A) result(A_inv) + real(rk), intent(in) :: A(:,:) + real(rk), allocatable :: A_inv(:,:) + + if (size(A,1) == size(A,2)) then + select case(size(A,1)) + case(2) + allocate(A_inv(size(A,1),size(A,2))) + A_inv(1,1) = A(2,2) + A_inv(1,2) = -A(1,2) + A_inv(2,1) = -A(2,1) + A_inv(2,2) = A(1,1) + A_inv = A_inv/det(A) + case(3) + allocate(A_inv(size(A,1),size(A,2))) + A_inv(1,1) = A(2,2)*A(3,3) - A(2,3)*A(3,2) + A_inv(1,2) = A(1,3)*A(3,2) - A(1,2)*A(3,3) + A_inv(1,3) = A(1,2)*A(2,3) - A(1,3)*A(2,2) + A_inv(2,1) = A(2,3)*A(3,1) - A(2,1)*A(3,3) + A_inv(2,2) = A(1,1)*A(3,3) - A(1,3)*A(3,1) + A_inv(2,3) = A(1,3)*A(2,1) - A(1,1)*A(2,3) + A_inv(3,1) = A(2,1)*A(3,2) - A(2,2)*A(3,1) + A_inv(3,2) = A(1,2)*A(3,1) - A(1,1)*A(3,2) + A_inv(3,3) = A(1,1)*A(2,2) - A(1,2)*A(2,1) + A_inv = A_inv/det(A) + end select + elseif (size(A,1)>size(A,2)) then + allocate(A_inv(size(A,2),size(A,1))) + A_inv = transpose(A) + A_inv = matmul(inv(matmul(A_inv, A)), A_inv) + elseif (size(A,1) author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function dyad_t1_t1(a, b) result(c) + real(rk), intent(in), contiguous :: a(:) + real(rk), intent(in), contiguous :: b(:) + real(rk), allocatable :: c(:,:) + integer :: i + + allocate(c(size(a), size(b))) + do concurrent(i = 1:size(c,1)) + c(i, :) = a(i) * b(:) + end do + end function + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine gauss_legendre_1D(interval, degree, Xksi, Wksi) + real(rk), intent(in) :: interval(2) + integer, intent(in) :: degree + real(rk), allocatable, intent(out) :: Xksi(:), Wksi(:) + + allocate(Xksi(degree+1), Wksi(degree+1)) + + call gauss_legendre(Xksi, Wksi, interval) + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine gauss_legendre_2D(interval1, interval2, degree, Xksi, Wksi) + real(rk), intent(in) :: interval1(2), interval2(2) + integer, intent(in) :: degree(2) + real(rk), allocatable, intent(out) :: Xksi(:,:), Wksi(:) + real(rk), allocatable :: Xksi1(:), Wksi1(:), Xksi2(:), Wksi2(:) + + allocate(Xksi1(degree(1)+1), Wksi1(degree(1)+1)) + allocate(Xksi2(degree(2)+1), Wksi2(degree(2)+1)) + + call gauss_legendre(Xksi1, Wksi1, interval1) + call gauss_legendre(Xksi2, Wksi2, interval2) + + call ndgrid(Xksi1, Xksi2, Xksi) + Wksi = kron(Wksi1, Wksi2) + end subroutine + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure subroutine gauss_legendre_3D(interval1, interval2, interval3, degree, Xksi, Wksi) + real(rk), intent(in) :: interval1(2), interval2(2), interval3(2) + integer, intent(in) :: degree(3) + real(rk), allocatable, intent(out) :: Xksi(:,:), Wksi(:) + real(rk), allocatable :: Xksi1(:), Wksi1(:), Xksi2(:), Wksi2(:), Xksi3(:), Wksi3(:) + + allocate(Xksi1(degree(1)+1), Wksi1(degree(1)+1)) + allocate(Xksi2(degree(2)+1), Wksi2(degree(2)+1)) + allocate(Xksi3(degree(3)+1), Wksi3(degree(3)+1)) + + call gauss_legendre(Xksi1, Wksi1, interval1) + call gauss_legendre(Xksi2, Wksi2, interval2) + call gauss_legendre(Xksi3, Wksi3, interval3) + + call ndgrid(Xksi1, Xksi2, Xksi3, Xksi) + Wksi = kron(kron(Wksi3, Wksi2), Wksi1) + end subroutine + !=============================================================================== + end module forcad_utils