Added cmp_elemFace*() and cmp_degreeFace() methods

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
 

ROADMAP.md

@@ -41,7 +41,8 @@ This roadmap outlines upcoming features and enhancements for ForCAD. Contributio
     - [ ] Improve PyVista scripts.
     - [ ] Use ForCAD for IGA simulations.
     - [ ] Use ForCAD for FEM simulations.
-    - [ ] Get surfaces and edges from volumes.
+    - [x] Get surfaces from volumes.
+    - [ ] Get edges from volumes.
     - [ ] Get edges from surfaces.
     - [ ] Add extrude method for surfaces.
     - [ ] Improve documentation.

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
     !-----------------------------------------------------------------------------

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,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
 
 !===============================================================================