From bb970073cf1c46bed9955c67d9e88b9ef72aad19 Mon Sep 17 00:00:00 2001 From: Seyed Ali Ghasemi Date: Tue, 21 May 2024 15:46:11 +0200 Subject: [PATCH] Added cmp_elemFace*() and cmp_degreeFace() methods - Update example_volume_1.f90. - Update CHANGELOG.md. --- CHANGELOG.md | 5 + example/example_volume_1.f90 | 23 ++++ src/forcad_nurbs_volume.f90 | 197 ++++++++++++++++++++++++++++++++++- 3 files changed, 223 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index df96a4139..3915ea6fc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -23,6 +23,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), - Added `Xt` to the `nurbs_surface` and `nurbs_volume` derived types. - Added codecoverage to the CI. - Added allocate(Tgci) and allocate(dTgci) to `compute_Tgc_nurbs_*d()` and `compute_dTgc_nurbs_*d()` subroutines. +- Added `cmp_elemFace_Xc_vis()` to extract element connectivity of the faces of control geometry. +- Added `cmp_elemFace_Xg_vis()` to extract element connectivity of the faces of visualization geometry points. +- Added `cmp_elemFace()` to extract element connectivity of faces. +- Added `cmp_degreeFace()` to extract the degrees of faces. +- Updated `example_volume_1` to use `cmp_elemFace()` and `cmp_degreeFace()`. ### Changed diff --git a/example/example_volume_1.f90 b/example/example_volume_1.f90 index 39fef0afc..920ce112f 100644 --- a/example/example_volume_1.f90 +++ b/example/example_volume_1.f90 @@ -145,6 +145,29 @@ program example3_volume !> Show the control geometry and geometry using PyVista call nurbs%show('vtk/nurbs_volume_Xc3.vtk','vtk/nurbs_volume_Xg3.vtk') + !----------------------------------------------------------------------------- + ! Extract faces + !----------------------------------------------------------------------------- + + !> first compute and set the connectivities of volume elements + call nurbs%set_elem(nurbs%cmp_elem()) + + !> get the connectivity of the face1 of the first element + print*, 'Face 1 of element 1:', nurbs%cmp_elemFace(elem=1, face=1) + print*, 'Face 2 of element 1:', nurbs%cmp_elemFace(elem=1, face=2) + print*, 'Face 3 of element 1:', nurbs%cmp_elemFace(elem=1, face=3) + print*, 'Face 4 of element 1:', nurbs%cmp_elemFace(elem=1, face=4) + print*, 'Face 5 of element 1:', nurbs%cmp_elemFace(elem=1, face=5) + print*, 'Face 6 of element 1:', nurbs%cmp_elemFace(elem=1, face=6) + + !> get the degree of the faces + print*, 'Degree of face 1:', nurbs%cmp_degreeFace(face=1) + print*, 'Degree of face 2:', nurbs%cmp_degreeFace(face=2) + print*, 'Degree of face 3:', nurbs%cmp_degreeFace(face=3) + print*, 'Degree of face 4:', nurbs%cmp_degreeFace(face=4) + print*, 'Degree of face 5:', nurbs%cmp_degreeFace(face=5) + print*, 'Degree of face 6:', nurbs%cmp_degreeFace(face=6) + !----------------------------------------------------------------------------- ! Finalizing !----------------------------------------------------------------------------- diff --git a/src/forcad_nurbs_volume.f90 b/src/forcad_nurbs_volume.f90 index cc9215ca7..174b2bb47 100644 --- a/src/forcad_nurbs_volume.f90 +++ b/src/forcad_nurbs_volume.f90 @@ -89,6 +89,12 @@ module forcad_nurbs_volume procedure :: show !!> Show the NURBS object using PyVista procedure :: nearest_point !!> Find the nearest point on the NURBS volume + ! Faces + procedure :: cmp_elemFace_Xc_vis !!> Compute faces of the control points + procedure :: cmp_elemFace_Xg_vis !!> Compute faces of the geometry points + procedure :: cmp_elemFace !!> Compute faces of the IGA elements + procedure :: cmp_degreeFace !!> Compute degrees of the faces + ! Shapes procedure :: set_hexahedron !!> Set a hexahedron procedure :: set_ring !!> Set a ring @@ -296,7 +302,7 @@ pure function compute_Xg_bspline_3d(f_Xt, f_knot1, f_knot2, f_knot3, f_degree, f if (allocated(this%Xg)) deallocate(this%Xg) allocate(this%Xg(this%ng(1)*this%ng(2)*this%ng(3), size(this%Xc,2))) - + if (allocated(this%Wc)) then ! NURBS volume this%Xg = compute_Xg_nurbs_3d(& this%Xt, this%knot1, this%knot2, this%knot3, this%degree, this%nc, this%ng, this%Xc, this%Wc) @@ -2602,7 +2608,7 @@ pure function nearest_point_help_3d(f_ng, f_Xg, f_point_Xg) result(f_distances) allocate(distances(this%ng(1)*this%ng(2)*this%ng(3))) distances = nearest_point_help_3d(this%ng, this%Xg, point_Xg) - + id_ = minloc(distances, dim=1) if (present(id)) id = id_ if (present(nearest_Xg)) nearest_Xg = this%Xg(id_,:) @@ -2610,6 +2616,193 @@ pure function nearest_point_help_3d(f_ng, f_Xg, f_point_Xg) result(f_distances) end subroutine !=============================================================================== + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function cmp_elemFace(this, elem, face) result(elemConn) + class(nurbs_volume), intent(in) :: this + integer, intent(in) :: elem + integer, intent(in) :: face + integer, allocatable :: elemConn(:) + integer :: n(3), ii, jj, k + + !> number of nodes in each direction + n = this%degree + 1 + + select case(face) + case(1) + allocate(elemConn(n(1)*n(2))) + elemConn = this%elemConn(elem, 1 : n(1)*n(2)) + case(2) + allocate(elemConn(n(1)*n(2))) + elemConn = this%elemConn(elem, n(1)*n(2)*n(3)-n(1)*n(2)+1 : n(1)*n(2)*n(3)) + case(3) + allocate(elemConn(n(1)*n(3))) + k = 0 + do jj = 1, n(3) + do ii = 1, n(1) + k = k+1 + elemConn(k) = this%elemConn(elem, n(1)*n(2)*jj + ii - n(1)*n(2)) + end do + end do + case(4) + allocate(elemConn(n(1)*n(3))) + k = 0 + do jj = 1, n(3) + do ii = 1, n(1) + k = k+1 + elemConn(k) = this%elemConn(elem, n(1)*n(2)*jj + ii - n(1)) + end do + end do + case(5) + allocate(elemConn(n(2)*n(3))) + do ii = 1, n(2)*n(3) + elemConn(ii) = this%elemConn(elem, n(1)*ii - n(1) + 1) + end do + case(6) + allocate(elemConn(n(2)*n(3))) + do ii = 1, n(2)*n(3) + elemConn(ii) = this%elemConn(elem, n(1)*ii) + end do + end select + end function + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function cmp_elemFace_Xc_vis(this, elem, face) result(elemConn) + class(nurbs_volume), intent(in) :: this + integer, intent(in) :: elem + integer, intent(in) :: face + integer, allocatable :: elemConn(:) + integer :: n(3), ii, jj, k + + !> number of nodes in each direction + n = [2,2,2] + + select case(face) + case(1) + allocate(elemConn(n(1)*n(2))) + elemConn = this%elemConn_Xc_vis(elem, 1 : n(1)*n(2)) + case(2) + allocate(elemConn(n(1)*n(2))) + elemConn = this%elemConn_Xc_vis(elem, n(1)*n(2)*n(3)-n(1)*n(2)+1 : n(1)*n(2)*n(3)) + case(3) + allocate(elemConn(n(1)*n(3))) + k = 0 + do jj = 1, n(3) + do ii = 1, n(1) + k = k+1 + elemConn(k) = this%elemConn_Xc_vis(elem, n(1)*n(2)*jj + ii - n(1)*n(2)) + end do + end do + case(4) + allocate(elemConn(n(1)*n(3))) + k = 0 + do jj = 1, n(3) + do ii = 1, n(1) + k = k+1 + elemConn(k) = this%elemConn_Xc_vis(elem, n(1)*n(2)*jj + ii - n(1)) + end do + end do + case(5) + allocate(elemConn(n(2)*n(3))) + do ii = 1, n(2)*n(3) + elemConn(ii) = this%elemConn_Xc_vis(elem, n(1)*ii - n(1) + 1) + end do + case(6) + allocate(elemConn(n(2)*n(3))) + do ii = 1, n(2)*n(3) + elemConn(ii) = this%elemConn_Xc_vis(elem, n(1)*ii) + end do + end select + end function + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function cmp_elemFace_Xg_vis(this, elem, face) result(elemConn) + class(nurbs_volume), intent(in) :: this + integer, intent(in) :: elem + integer, intent(in) :: face + integer, allocatable :: elemConn(:) + integer :: n(3), ii, jj, k + + !> number of nodes in each direction + n = [2,2,2] + + select case(face) + case(1) + allocate(elemConn(n(1)*n(2))) + elemConn = this%elemConn_Xg_vis(elem, 1 : n(1)*n(2)) + case(2) + allocate(elemConn(n(1)*n(2))) + elemConn = this%elemConn_Xg_vis(elem, n(1)*n(2)*n(3)-n(1)*n(2)+1 : n(1)*n(2)*n(3)) + case(3) + allocate(elemConn(n(1)*n(3))) + k = 0 + do jj = 1, n(3) + do ii = 1, n(1) + k = k+1 + elemConn(k) = this%elemConn_Xg_vis(elem, n(1)*n(2)*jj + ii - n(1)*n(2)) + end do + end do + case(4) + allocate(elemConn(n(1)*n(3))) + k = 0 + do jj = 1, n(3) + do ii = 1, n(1) + k = k+1 + elemConn(k) = this%elemConn_Xg_vis(elem, n(1)*n(2)*jj + ii - n(1)) + end do + end do + case(5) + allocate(elemConn(n(2)*n(3))) + do ii = 1, n(2)*n(3) + elemConn(ii) = this%elemConn_Xg_vis(elem, n(1)*ii - n(1) + 1) + end do + case(6) + allocate(elemConn(n(2)*n(3))) + do ii = 1, n(2)*n(3) + elemConn(ii) = this%elemConn_Xg_vis(elem, n(1)*ii) + end do + end select + end function + !=============================================================================== + + + !=============================================================================== + !> author: Seyed Ali Ghasemi + !> license: BSD 3-Clause + pure function cmp_degreeFace(this, face) result(degree) + class(nurbs_volume), intent(in) :: this + integer, intent(in) :: face + integer :: degree(3) + + select case (face) + case(1) + degree = [this%degree(1), this%degree(2), 0] + case(2) + degree = [this%degree(1), this%degree(2), 0] + case(3) + degree = [this%degree(1), 0, this%degree(3)] + case(4) + degree = [this%degree(1), 0, this%degree(3)] + case(5) + degree = [0, this%degree(2), this%degree(3)] + case(6) + degree = [0, this%degree(2), this%degree(3)] + case default + error stop 'Invalid face number' + end select + end function + !=============================================================================== + end module forcad_nurbs_volume !===============================================================================