make_annulus Subroutine

public subroutine make_annulus(mesh, radius, inner_radius, n_theta, n_r, center)

XY平面上の同心リングを極座標分割し、三角形メッシュを生成する。

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 :: inner_radius

外半径 [m]（省略時 0.5）。 内半径 [m]（省略時 0.25）。
integer(kind=i32), intent(in), optional :: n_theta

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

周方向分割数（省略時 24）。 半径方向分割数（省略時 4）。
real(kind=dp), intent(in), optional :: center(3)

外半径 [m]（省略時 0.5）。 内半径 [m]（省略時 0.25）。

Calls

proc~~make_annulus~~CallsGraph proc~make_annulus make_annulus proc~init_mesh init_mesh proc~make_annulus->proc~init_mesh proc~push_tri push_tri proc~make_annulus->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_annulus~~CalledByGraph proc~make_annulus make_annulus proc~build_one_template build_one_template proc~build_one_template->proc~make_annulus proc~make_disk make_disk proc~build_one_template->proc~make_disk proc~make_disk->proc~make_annulus 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_annulus(mesh, radius, inner_radius, n_theta, n_r, center)
    type(mesh_type), intent(out) :: mesh
    real(dp), intent(in), optional :: radius, inner_radius, center(3)
    integer(i32), intent(in), optional :: n_theta, n_r
    real(dp) :: r_out, r_in, c(3), pi, t0, t1, r0, r1, dr
    integer(i32) :: nt, nr0, ir, it, itri, nelem
    real(dp) :: p00(3), p01(3), p10(3), p11(3), center_p(3)
    real(dp), allocatable :: v0(:, :), v1(:, :), v2(:, :)

    r_out = 0.5d0; r_in = 0.25d0; nt = 24; nr0 = 4; c = 0.0d0
    if (present(radius)) r_out = radius
    if (present(inner_radius)) r_in = inner_radius
    if (present(n_theta)) nt = n_theta
    if (present(n_r)) nr0 = n_r
    if (present(center)) c = center
    if (nt < 3 .or. nr0 <= 0) error stop "n_theta >= 3 and n_r > 0"
    if (r_out <= 0.0d0) error stop "radius must be > 0"
    if (r_in < 0.0d0) error stop "inner_radius must be >= 0"
    if (r_in >= r_out) error stop "inner_radius must be smaller than radius"

    if (r_in <= 0.0d0) then
      nelem = nt*(2*nr0 - 1)
    else
      nelem = 2*nt*nr0
    end if
    allocate (v0(3, nelem), v1(3, nelem), v2(3, nelem))

    pi = acos(-1.0d0)
    dr = (r_out - r_in)/real(nr0, dp)
    center_p = [c(1), c(2), c(3)]
    itri = 0
    do ir = 0, nr0 - 1
      r0 = r_in + dr*real(ir, dp)
      r1 = r_in + dr*real(ir + 1, 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)
        p00 = [c(1) + r0*cos(t0), c(2) + r0*sin(t0), c(3)]
        p01 = [c(1) + r0*cos(t1), c(2) + r0*sin(t1), c(3)]
        p10 = [c(1) + r1*cos(t0), c(2) + r1*sin(t0), c(3)]
        p11 = [c(1) + r1*cos(t1), c(2) + r1*sin(t1), c(3)]
        if (ir == 0 .and. r_in <= 0.0d0) then
          call push_tri(v0, v1, v2, itri, center_p, p10, p11)
        else
          call push_tri(v0, v1, v2, itri, p00, p10, p11)
          call push_tri(v0, v1, v2, itri, p00, p11, p01)
        end if
      end do
    end do

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