make_cylinder Subroutine

public subroutine make_cylinder(mesh, radius, height, n_theta, n_z, cap, center, cap_top, cap_bottom, uniform_cap, n_r)

円柱側面を分割生成し、必要に応じて上下面キャップを追加したメッシュを生成する。 （省略時 .false.: 中心から外周への n_theta 枚の扇）。

Arguments

Type IntentOptional Attributes Name
type(mesh_type), intent(out) :: mesh

生成した円柱メッシュ。
real(kind=dp), intent(in), optional :: radius

円柱半径 [m]（省略時 0.5）。
real(kind=dp), intent(in), optional :: height

円柱半径 [m]（省略時 0.5）。 円柱高さ [m]（省略時 1.0）。
integer(kind=i32), intent(in), optional :: n_theta

周方向分割数（省略時 24）。
integer(kind=i32), intent(in), optional :: n_z

周方向分割数（省略時 24）。 軸方向分割数（省略時 1）。
logical, intent(in), optional :: cap

上下面キャップをまとめて指定する後方互換フラグ（省略時 .true.）。
real(kind=dp), intent(in), optional :: center(3)

円柱半径 [m]（省略時 0.5）。 円柱高さ [m]（省略時 1.0）。
logical, intent(in), optional :: cap_top

上下面キャップをまとめて指定する後方互換フラグ（省略時 .true.）。 上面キャップを生成するか（省略時 .true.）。
logical, intent(in), optional :: cap_bottom

上下面キャップをまとめて指定する後方互換フラグ（省略時 .true.）。 上面キャップを生成するか（省略時 .true.）。 下面キャップを生成するか（省略時 .true.）。
logical, intent(in), optional :: uniform_cap

上下面キャップをまとめて指定する後方互換フラグ（省略時 .true.）。 上面キャップを生成するか（省略時 .true.）。 下面キャップを生成するか（省略時 .true.）。 .true. ならキャップを n_r リングの一様配置の円板にする
integer(kind=i32), intent(in), optional :: n_r

周方向分割数（省略時 24）。 軸方向分割数（省略時 1）。 uniform_cap のキャップの半径方向分割数（省略時 4）。

Calls

proc~~make_cylinder~~CallsGraph proc~make_cylinder make_cylinder proc~init_mesh init_mesh proc~make_cylinder->proc~init_mesh proc~push_tri push_tri proc~make_cylinder->proc~push_tri proc~update_mesh_geometry update_mesh_geometry proc~init_mesh->proc~update_mesh_geometry proc~build_collision_grid build_collision_grid proc~update_mesh_geometry->proc~build_collision_grid proc~cross~2 cross proc~update_mesh_geometry->proc~cross~2 proc~fill_panel_quadrature fill_panel_quadrature proc~update_mesh_geometry->proc~fill_panel_quadrature proc~init_panel_geometry init_panel_geometry proc~update_mesh_geometry->proc~init_panel_geometry proc~cell_id cell_id proc~build_collision_grid->proc~cell_id proc~coord_to_cell coord_to_cell proc~build_collision_grid->proc~coord_to_cell

Called by

proc~~make_cylinder~~CalledByGraph proc~make_cylinder make_cylinder proc~build_one_template build_one_template proc~build_one_template->proc~make_cylinder proc~build_parts_mesh build_parts_mesh proc~build_parts_mesh->proc~build_one_template proc~build_mesh_from_config build_mesh_from_config proc~build_mesh_from_config->proc~build_parts_mesh proc~load_or_init_run_state load_or_init_run_state proc~load_or_init_run_state->proc~build_mesh_from_config program~main main program~main->proc~load_or_init_run_state

Source Code

  subroutine make_cylinder(mesh, radius, height, n_theta, n_z, cap, center, cap_top, cap_bottom, uniform_cap, n_r)
    type(mesh_type), intent(out) :: mesh
    real(dp), intent(in), optional :: radius, height, center(3)
    integer(i32), intent(in), optional :: n_theta, n_z, n_r
    logical, intent(in), optional :: cap, cap_top, cap_bottom, uniform_cap
    real(dp) :: r, h, c(3), z0, z1
    integer(i32) :: nt, nz0, iz, it, nelem, itri, ncap, nr0
    logical :: cap_top0, cap_bottom0, uniform_cap0
    real(dp), allocatable :: v0(:, :), v1(:, :), v2(:, :), circle(:, :)

    r = 0.5d0; h = 1.0d0; nt = 24; nz0 = 1; cap_top0 = .true.; cap_bottom0 = .true.; c = 0.0d0
    if (present(radius)) r = radius
    if (present(height)) h = height
    if (present(n_theta)) nt = n_theta
    if (present(n_z)) nz0 = n_z
    if (present(cap)) then
      cap_top0 = cap
      cap_bottom0 = cap
    end if
    if (present(cap_top)) cap_top0 = cap_top
    if (present(cap_bottom)) cap_bottom0 = cap_bottom
    if (present(center)) c = center
    uniform_cap0 = .false.; nr0 = 4
    if (present(uniform_cap)) uniform_cap0 = uniform_cap
    if (present(n_r)) nr0 = n_r
    if (nt < 3 .or. nz0 <= 0) error stop "n_theta >= 3 and n_z > 0"
    if (uniform_cap0 .and. nr0 <= 0) error stop "n_r > 0"

    ncap = nt
    if (uniform_cap0) ncap = uniform_annulus_count(0.0d0, r, nt, nr0)
    nelem = 2*nt*nz0
    if (cap_bottom0) nelem = nelem + ncap
    if (cap_top0) nelem = nelem + ncap
    allocate (v0(3, nelem), v1(3, nelem), v2(3, nelem)); itri = 0

    call build_unit_circle(nt, circle)
    do iz = 0, nz0 - 1
      z0 = c(3) - 0.5d0*h + h*real(iz, dp)/real(nz0, dp)
      z1 = c(3) - 0.5d0*h + h*real(iz + 1, dp)/real(nz0, dp)
      ! The top ring uses the cap height expression so side and top cap share their rim vertices.
      if (iz == nz0 - 1) z1 = c(3) + 0.5d0*h
      do it = 0, nt - 1
        call push_tri(v0, v1, v2, itri, &
                      [c(1) + r*circle(1, it), c(2) + r*circle(2, it), z0], &
                      [c(1) + r*circle(1, it + 1), c(2) + r*circle(2, it + 1), z0], &
                      [c(1) + r*circle(1, it + 1), c(2) + r*circle(2, it + 1), z1])
        call push_tri(v0, v1, v2, itri, &
                      [c(1) + r*circle(1, it), c(2) + r*circle(2, it), z0], &
                      [c(1) + r*circle(1, it + 1), c(2) + r*circle(2, it + 1), z1], &
                      [c(1) + r*circle(1, it), c(2) + r*circle(2, it), z1])
      end do
    end do

    ! A uniform cap's outermost ring reuses the side circle, so the cap rim sits on the side vertices.
    if (cap_bottom0) then
      if (uniform_cap0) then
        call push_uniform_annulus(v0, v1, v2, itri, [c(1), c(2), c(3) - 0.5d0*h], 0.0d0, r, nt, nr0, .true.)
      else
        do it = 0, nt - 1
          call push_tri(v0, v1, v2, itri, &
                        [c(1), c(2), c(3) - 0.5d0*h], &
                        [c(1) + r*circle(1, it + 1), c(2) + r*circle(2, it + 1), c(3) - 0.5d0*h], &
                        [c(1) + r*circle(1, it), c(2) + r*circle(2, it), c(3) - 0.5d0*h])
        end do
      end if
    end if
    if (cap_top0) then
      if (uniform_cap0) then
        call push_uniform_annulus(v0, v1, v2, itri, [c(1), c(2), c(3) + 0.5d0*h], 0.0d0, r, nt, nr0, .false.)
      else
        do it = 0, nt - 1
          call push_tri(v0, v1, v2, itri, &
                        [c(1), c(2), c(3) + 0.5d0*h], &
                        [c(1) + r*circle(1, it), c(2) + r*circle(2, it), c(3) + 0.5d0*h], &
                        [c(1) + r*circle(1, it + 1), c(2) + r*circle(2, it + 1), c(3) + 0.5d0*h])
        end do
      end if
    end if

    call init_mesh(mesh, v0, v1, v2)
  end subroutine make_cylinder