Format and analytical levelset refactor

haochey · haochey · commit 33c314aae1b1 · 2024-11-18T17:35:53.000-05:00
diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp
@@ -35,6 +35,6 @@ module m_constants
     real(kind(0d0)), parameter :: threshold_vector_zero = 1d-10 !< Threshold to treat the component of a vector to be zero
     real(kind(0d0)), parameter :: threshold_edge_zero = 1d-8 !< Threshold to treat two edges to be overlapped
     real(kind(0d0)), parameter :: threshold_bary = 1d-1 !< Threshold to interpolate a barycentric facet
-    real(kind(0d0)), parameter :: initial_distance_buffer = 1d12 !< Initialized levelset distance for the shortest path pair algorithm 
+    real(kind(0d0)), parameter :: initial_distance_buffer = 1d12 !< Initialized levelset distance for the shortest path pair algorithm
 
 end module m_constants
diff --git a/src/pre_process/m_compute_levelset.fpp b/src/pre_process/m_compute_levelset.fpp
@@ -20,12 +20,12 @@ module m_compute_levelset
 
     implicit none
 
-    private; public :: s_compute_cylinder_levelset, s_compute_circle_levelset, &
- s_compute_airfoil_levelset, &
- s_compute_3D_airfoil_levelset, &
- s_compute_rectangle_levelset, &
- s_compute_cuboid_levelset, &
- s_compute_sphere_levelset
+    private; public :: s_cylinder_levelset, s_circle_levelset, &
+ s_airfoil_levelset, &
+ s_3D_airfoil_levelset, &
+ s_rectangle_levelset, &
+ s_cuboid_levelset, &
+ s_sphere_levelset
 
     real(kind(0d0)) :: x_centroid, y_centroid, z_centroid
     real(kind(0d0)) :: length_x, length_y, length_z
@@ -39,7 +39,7 @@ module m_compute_levelset
 
 contains
 
-    subroutine s_compute_circle_levelset(levelset, levelset_norm, ib_patch_id)
+    subroutine s_circle_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -73,9 +73,9 @@ contains
             end do
         end do
 
-    end subroutine s_compute_circle_levelset
+    end subroutine s_circle_levelset
 
-    subroutine s_compute_airfoil_levelset(levelset, levelset_norm, ib_patch_id)
+    subroutine s_airfoil_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -156,9 +156,9 @@ contains
             end do
         end do
 
-    end subroutine s_compute_airfoil_levelset
+    end subroutine s_airfoil_levelset
 
-    subroutine s_compute_3D_airfoil_levelset(levelset, levelset_norm, ib_patch_id)
+    subroutine s_3D_airfoil_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -257,10 +257,10 @@ contains
             end do
         end do
 
-    end subroutine s_compute_3D_airfoil_levelset
+    end subroutine s_3D_airfoil_levelset
 
     !>  Initialize IBM module
-    subroutine s_compute_rectangle_levelset(levelset, levelset_norm, ib_patch_id)
+    subroutine s_rectangle_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -345,9 +345,9 @@ contains
             end do
         end do
 
-    end subroutine s_compute_rectangle_levelset
+    end subroutine s_rectangle_levelset
 
-    subroutine s_compute_cuboid_levelset(levelset, levelset_norm, ib_patch_id)
+    subroutine s_cuboid_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -459,9 +459,9 @@ contains
             end do
         end do
 
-    end subroutine s_compute_cuboid_levelset
+    end subroutine s_cuboid_levelset
 
-    subroutine s_compute_sphere_levelset(levelset, levelset_norm, ib_patch_id)
+    subroutine s_sphere_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -496,9 +496,9 @@ contains
             end do
         end do
 
-    end subroutine s_compute_sphere_levelset
+    end subroutine s_sphere_levelset
 
-    subroutine s_compute_cylinder_levelset(levelset, levelset_norm, ib_patch_id)
+    subroutine s_cylinder_levelset(levelset, levelset_norm, ib_patch_id)
 
         type(levelset_field), intent(INOUT) :: levelset
         type(levelset_norm_field), intent(INOUT) :: levelset_norm
@@ -569,6 +569,6 @@ contains
             end do
         end do
 
-    end subroutine s_compute_cylinder_levelset
+    end subroutine s_cylinder_levelset
 
 end module m_compute_levelset
diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp
@@ -149,15 +149,15 @@ contains
                 !> @{
                 ! Spherical patch
                 if (patch_icpp(i)%geometry == 8) then
-                    call s_sphere(i, patch_id_fp, q_prim_vf, .false.)
+                    call s_sphere(i, patch_id_fp, q_prim_vf)
 
                     ! Cuboidal patch
                 elseif (patch_icpp(i)%geometry == 9) then
                     call s_cuboid(i, patch_id_fp, q_prim_vf)
 
                     ! Cylindrical patch
                 elseif (patch_icpp(i)%geometry == 10) then
-                    call s_cylinder(i, patch_id_fp, q_prim_vf, .false.)
+                    call s_cylinder(i, patch_id_fp, q_prim_vf)
 
                     ! Swept plane patch
                 elseif (patch_icpp(i)%geometry == 11) then
@@ -181,7 +181,7 @@ contains
 
                     ! 3D STL patch
                 elseif (patch_icpp(i)%geometry == 21) then
-                    call s_model(i, patch_id_fp, q_prim_vf, .false.)
+                    call s_model(i, patch_id_fp, q_prim_vf)
 
                 end if
 
@@ -197,19 +197,19 @@ contains
                 end if
 
                 if (patch_ib(i)%geometry == 8) then
-                    call s_sphere(i, ib_markers%sf, q_prim_vf, .true.)
-                    call s_compute_sphere_levelset(levelset, levelset_norm, i)
+                    call s_sphere(i, ib_markers%sf, q_prim_vf, ib)
+                    call s_sphere_levelset(levelset, levelset_norm, i)
                     ! Cylindrical patch
                 elseif (patch_ib(i)%geometry == 10) then
-                    call s_cylinder(i, ib_markers%sf, q_prim_vf, .true.)
-                    call s_compute_cylinder_levelset(levelset, levelset_norm, i)
+                    call s_cylinder(i, ib_markers%sf, q_prim_vf, ib)
+                    call s_cylinder_levelset(levelset, levelset_norm, i)
                 elseif (patch_ib(i)%geometry == 11) then
-                    call s_3D_airfoil(i, ib_markers%sf, q_prim_vf, .true.)
-                    call s_compute_3D_airfoil_levelset(levelset, levelset_norm, i)
+                    call s_3D_airfoil(i, ib_markers%sf, q_prim_vf, ib)
+                    call s_3D_airfoil_levelset(levelset, levelset_norm, i)
 
                     ! STL+IBM patch
                 elseif (patch_ib(i)%geometry == 12) then
-                    call s_model(i, ib_markers%sf, q_prim_vf, .true., levelset, levelset_norm)
+                    call s_model(i, ib_markers%sf, q_prim_vf, ib, levelset, levelset_norm)
                 end if
             end do
             !> @}
@@ -279,18 +279,18 @@ contains
                     print *, 'Processing 2D ib patch ', i
                 end if
                 if (patch_ib(i)%geometry == 2) then
-                    call s_circle(i, ib_markers%sf, q_prim_vf, .true.)
-                    call s_compute_circle_levelset(levelset, levelset_norm, i)
+                    call s_circle(i, ib_markers%sf, q_prim_vf, ib)
+                    call s_circle_levelset(levelset, levelset_norm, i)
                     ! Rectangular patch
                 elseif (patch_ib(i)%geometry == 3) then
-                    call s_rectangle(i, ib_markers%sf, q_prim_vf, .true.)
-                    call s_compute_rectangle_levelset(levelset, levelset_norm, i)
+                    call s_rectangle(i, ib_markers%sf, q_prim_vf, ib)
+                    call s_rectangle_levelset(levelset, levelset_norm, i)
                 elseif (patch_ib(i)%geometry == 4) then
-                    call s_airfoil(i, ib_markers%sf, q_prim_vf, .true.)
-                    call s_compute_airfoil_levelset(levelset, levelset_norm, i)
+                    call s_airfoil(i, ib_markers%sf, q_prim_vf, ib)
+                    call s_airfoil_levelset(levelset, levelset_norm, i)
                     ! STL+IBM patch
                 elseif (patch_ib(i)%geometry == 5) then
-                    call s_model(i, ib_markers%sf, q_prim_vf, .true., levelset, levelset_norm)
+                    call s_model(i, ib_markers%sf, q_prim_vf, ib, levelset, levelset_norm)
                 end if
             end do
             !> @}
diff --git a/src/pre_process/m_model.fpp b/src/pre_process/m_model.fpp
@@ -704,13 +704,13 @@ contains
         interpolate = .false.
 
         do j = 1, boundary_edge_count
-            edge_v(1,1) = boundary_v(j, 1, 1)
-            edge_v(2,1) = boundary_v(j, 2, 1)
-            edge_v(1,2) = boundary_v(j, 1, 2)
-            edge_v(2,2) = boundary_v(j, 2, 2)
+            edge_v(1, 1) = boundary_v(j, 1, 1)
+            edge_v(2, 1) = boundary_v(j, 2, 1)
+            edge_v(1, 2) = boundary_v(j, 1, 2)
+            edge_v(2, 2) = boundary_v(j, 2, 2)
 
-            l1 = dsqrt((edge_v(2,1) - edge_v(1,1))**2 + &
-                       (edge_v(2,2) - edge_v(1,2))**2)
+            l1 = dsqrt((edge_v(2, 1) - edge_v(1, 1))**2 + &
+                       (edge_v(2, 2) - edge_v(1, 2))**2)
 
             if ((l1 > cell_width)) then
                 interpolate = .true.
@@ -737,23 +737,23 @@ contains
 
         do i = 1, model%ntrs
             do j = 1, 3
-               tri_v(1, j) = model%trs(i)%v(1, j)
-               tri_v(2, j) = model%trs(i)%v(2, j)
-               tri_v(3, j) = model%trs(i)%v(3, j)
+                tri_v(1, j) = model%trs(i)%v(1, j)
+                tri_v(2, j) = model%trs(i)%v(2, j)
+                tri_v(3, j) = model%trs(i)%v(3, j)
             end do
 
-            edge_l(1) = dsqrt((tri_v(1, 2) - tri_v(1, 1))**2 + & 
-                              (tri_v(2, 2) - tri_v(2, 1))**2 + & 
+            edge_l(1) = dsqrt((tri_v(1, 2) - tri_v(1, 1))**2 + &
+                              (tri_v(2, 2) - tri_v(2, 1))**2 + &
                               (tri_v(3, 2) - tri_v(3, 1))**2)
-            edge_l(2) = dsqrt((tri_v(1, 3) - tri_v(1, 2))**2 + & 
-                              (tri_v(2, 3) - tri_v(2, 2))**2 + & 
+            edge_l(2) = dsqrt((tri_v(1, 3) - tri_v(1, 2))**2 + &
+                              (tri_v(2, 3) - tri_v(2, 2))**2 + &
                               (tri_v(3, 3) - tri_v(3, 2))**2)
-            edge_l(3) = dsqrt((tri_v(1, 1) - tri_v(1, 3))**2 + & 
-                              (tri_v(2, 1) - tri_v(2, 3))**2 + & 
+            edge_l(3) = dsqrt((tri_v(1, 1) - tri_v(1, 3))**2 + &
+                              (tri_v(2, 1) - tri_v(2, 3))**2 + &
                               (tri_v(3, 1) - tri_v(3, 3))**2)
 
             if ((edge_l(1) > cell_width) .or. &
-                (edge_l(2) > cell_width) .or. & 
+                (edge_l(2) > cell_width) .or. &
                 (edge_l(3) > cell_width)) then
                 interpolate = .true.
             end if
@@ -791,7 +791,7 @@ contains
             edge_y(2) = boundary_v(i, 2, 2)
 
             ! Compute the length of the edge
-            edge_length = dsqrt((edge_x(2) - edge_x(1))**2 + & 
+            edge_length = dsqrt((edge_x(2) - edge_x(1))**2 + &
                                 (edge_y(2) - edge_y(1))**2)
 
             ! Determine the number of segments
@@ -818,7 +818,7 @@ contains
             edge_y(2) = boundary_v(i, 2, 2)
 
             ! Compute the length of the edge
-            edge_length = dsqrt((edge_x(2) - edge_x(1))**2 + & 
+            edge_length = dsqrt((edge_x(2) - edge_x(1))**2 + &
                                 (edge_y(2) - edge_y(1))**2)
 
             ! Determine the number of segments and interpolation step
@@ -896,7 +896,7 @@ contains
                 tri(2, 3) = model%trs(i)%v(mod(j, 3) + 1, 3)
 
                 ! Compute the length of the edge
-                edge_length = dsqrt((tri(2, 1) - tri(1, 1))**2 + & 
+                edge_length = dsqrt((tri(2, 1) - tri(1, 1))**2 + &
                                     (tri(2, 2) - tri(1, 2))**2 + &
                                     (tri(2, 3) - tri(1, 3))**2)
 
@@ -943,7 +943,7 @@ contains
                 tri(2, 3) = model%trs(i)%v(mod(j, 3) + 1, 3)
 
                 ! Compute the length of the edge
-                edge_length = dsqrt((tri(2, 1) - tri(1, 1))**2 + & 
+                edge_length = dsqrt((tri(2, 1) - tri(1, 1))**2 + &
                                     (tri(2, 2) - tri(1, 2))**2 + &
                                     (tri(2, 3) - tri(1, 3))**2)
 
@@ -1034,13 +1034,13 @@ contains
                 tri(j, 2) = model%trs(i)%v(j, 2)
                 tri(j, 3) = model%trs(i)%v(j, 3)
                 dist_buffer(j) = dsqrt((point(1) - tri(j, 1))**2 + &
-                                    (point(2) - tri(j, 2))**2 + &
-                                    (point(3) - tri(j, 3))**2)
+                                       (point(2) - tri(j, 2))**2 + &
+                                       (point(3) - tri(j, 3))**2)
             end do
-            
+
             ! Get the surface center of each triangle facet
             do j = 1, 3
-                midp(j) = (tri(1, j) + tri(2, j)+ tri(3, j))/3
+                midp(j) = (tri(1, j) + tri(2, j) + tri(3, j))/3
             end do
 
             dist_t_min = minval(dist_buffer(1:3))
@@ -1163,7 +1163,7 @@ contains
         do i = 1, total_vertices
             dist_buffer = dsqrt((point(1) - interpolated_boundary_v(i, 1))**2 + &
                                 (point(2) - interpolated_boundary_v(i, 2))**2 + &
-                                (point(3)- interpolated_boundary_v(i, 3))**2)
+                                (point(3) - interpolated_boundary_v(i, 3))**2)
 
             if (min_dist > dist_buffer) then
                 min_dist = dist_buffer
diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp
