make_cylinder Subroutine

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

円柱側面を分割生成し、必要に応じて上下面キャップを追加したメッシュを生成する。

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.）。

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~cell_id~2 cell_id proc~build_collision_grid->proc~cell_id~2 proc~coord_to_cell~2 coord_to_cell proc~build_collision_grid->proc~coord_to_cell~2

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_template_mesh build_template_mesh proc~build_template_mesh->proc~build_one_template proc~build_mesh_from_config build_mesh_from_config proc~build_mesh_from_config->proc~build_template_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)
    type(mesh_type), intent(out) :: mesh
    real(dp), intent(in), optional :: radius, height, center(3)
    integer(i32), intent(in), optional :: n_theta, n_z
    logical, intent(in), optional :: cap, cap_top, cap_bottom
    real(dp) :: r, h, c(3), t0, t1, z0, z1, pi
    integer(i32) :: nt, nz0, iz, it, nelem, itri
    logical :: cap_top0, cap_bottom0
    real(dp), allocatable :: v0(:, :), v1(:, :), v2(:, :)

    pi = acos(-1.0d0)
    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
    if (nt < 3 .or. nz0 <= 0) error stop "n_theta >= 3 and n_z > 0"

    nelem = 2*nt*nz0
    if (cap_bottom0) nelem = nelem + nt
    if (cap_top0) nelem = nelem + nt
    allocate (v0(3, nelem), v1(3, nelem), v2(3, nelem)); itri = 0

    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)
      do it = 0, nt - 1
        t0 = 2.0d0*pi*real(it, dp)/real(nt, dp)
        t1 = 2.0d0*pi*real(mod(it + 1, nt), dp)/real(nt, dp)
        call push_tri(v0, v1, v2, itri, &
                      [c(1) + r*cos(t0), c(2) + r*sin(t0), z0], &
                      [c(1) + r*cos(t1), c(2) + r*sin(t1), z0], &
                      [c(1) + r*cos(t1), c(2) + r*sin(t1), z1])
        call push_tri(v0, v1, v2, itri, &
                      [c(1) + r*cos(t0), c(2) + r*sin(t0), z0], &
                      [c(1) + r*cos(t1), c(2) + r*sin(t1), z1], &
                      [c(1) + r*cos(t0), c(2) + r*sin(t0), z1])
      end do
    end do

    if (cap_bottom0) then
      do it = 0, nt - 1
        t0 = 2.0d0*pi*real(it, dp)/real(nt, dp)
        t1 = 2.0d0*pi*real(mod(it + 1, nt), dp)/real(nt, dp)
        call push_tri(v0, v1, v2, itri, &
                      [c(1), c(2), c(3) - 0.5d0*h], &
                      [c(1) + r*cos(t1), c(2) + r*sin(t1), c(3) - 0.5d0*h], &
                      [c(1) + r*cos(t0), c(2) + r*sin(t0), c(3) - 0.5d0*h])
      end do
    end if
    if (cap_top0) then
      do it = 0, nt - 1
        t0 = 2.0d0*pi*real(it, dp)/real(nt, dp)
        t1 = 2.0d0*pi*real(mod(it + 1, nt), dp)/real(nt, dp)
        call push_tri(v0, v1, v2, itri, &
                      [c(1), c(2), c(3) + 0.5d0*h], &
                      [c(1) + r*cos(t0), c(2) + r*sin(t0), c(3) + 0.5d0*h], &
                      [c(1) + r*cos(t1), c(2) + r*sin(t1), c(3) + 0.5d0*h])
      end do
    end if

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