Skip to content

Commit

Permalink
Added cmp_elemFace*() and cmp_degreeFace() methods
Browse files Browse the repository at this point in the history 
- Update example_volume_1.f90.
- Update CHANGELOG.md.
  • Loading branch information
@gha3mi
gha3mi committed May 21, 2024
1 parent 9422e8e commit bb97007
Show file tree
Hide file tree
Showing 3 changed files with 223 additions and 2 deletions.
5 changes: 5 additions & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
23 changes: 23 additions & 0 deletions example/example_volume_1.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
!-----------------------------------------------------------------------------
Expand Down
197 changes: 195 additions & 2 deletions src/forcad_nurbs_volume.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -2602,14 +2608,201 @@ 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_,:)
if (present(nearest_Xt)) nearest_Xt = this%Xt(id_,:)
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

!===============================================================================
Expand Down

0 comments on commit bb97007

Please sign in to comment.