make_sphere Subroutine

public subroutine make_sphere(mesh, radius, n_lon, n_lat, center)

経度・緯度分割に基づき球面三角形メッシュを生成する。

Arguments

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

生成した球面三角形メッシュ。
real(kind=dp), intent(in), optional :: radius

球半径 [m]（省略時 0.5）。
integer(kind=i32), intent(in), optional :: n_lon

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

経度方向分割数（省略時 24）。 緯度方向分割数（省略時 12）。
real(kind=dp), intent(in), optional :: center(3)

球半径 [m]（省略時 0.5）。

Calls

proc~~make_sphere~~CallsGraph proc~make_sphere make_sphere proc~init_mesh init_mesh proc~make_sphere->proc~init_mesh proc~push_tri push_tri proc~make_sphere->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~3 cross proc~update_mesh_geometry->proc~cross~3 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~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_sphere~~CalledByGraph proc~make_sphere make_sphere proc~build_one_template build_one_template proc~build_one_template->proc~make_sphere 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_sphere(mesh, radius, n_lon, n_lat, center)
    type(mesh_type), intent(out) :: mesh
    real(dp), intent(in), optional :: radius, center(3)
    integer(i32), intent(in), optional :: n_lon, n_lat
    real(dp) :: r, c(3), p0, p1, pi, ring_r0, ring_r1, ring_z0, ring_z1
    integer(i32) :: nl, nphi, ilon, ilat, nelem, itri
    real(dp) :: a(3), b(3), c0(3), d(3)
    real(dp), allocatable :: v0(:, :), v1(:, :), v2(:, :), circle(:, :)

    pi = acos(-1.0d0)
    r = 0.5d0; nl = 24; nphi = 12; c = 0.0d0
    if (present(radius)) r = radius
    if (present(n_lon)) nl = n_lon
    if (present(n_lat)) nphi = n_lat
    if (present(center)) c = center
    if (nl < 3 .or. nphi < 2) error stop "n_lon >= 3 and n_lat >= 2"

    nelem = 2*nl*(nphi - 1)
    allocate (v0(3, nelem), v1(3, nelem), v2(3, nelem)); itri = 0

    call build_unit_circle(nl, circle)
    do ilat = 0, nphi - 1
      p0 = pi*real(ilat, dp)/real(nphi, dp)
      p1 = pi*real(ilat + 1, dp)/real(nphi, dp)
      ring_r0 = r*sin(p0)
      ring_r1 = r*sin(p1)
      ring_z0 = c(3) + r*cos(p0)
      ring_z1 = c(3) + r*cos(p1)
      do ilon = 0, nl - 1
        a = [c(1) + ring_r0*circle(1, ilon), c(2) + ring_r0*circle(2, ilon), ring_z0]
        b = [c(1) + ring_r0*circle(1, ilon + 1), c(2) + ring_r0*circle(2, ilon + 1), ring_z0]
        c0 = [c(1) + ring_r1*circle(1, ilon), c(2) + ring_r1*circle(2, ilon), ring_z1]
        d = [c(1) + ring_r1*circle(1, ilon + 1), c(2) + ring_r1*circle(2, ilon + 1), ring_z1]
        if (ilat == 0) then
          call push_tri(v0, v1, v2, itri, a, c0, d)
        else if (ilat == nphi - 1) then
          call push_tri(v0, v1, v2, itri, a, d, b)
        else
          call push_tri(v0, v1, v2, itri, a, d, b)
          call push_tri(v0, v1, v2, itri, a, c0, d)
        end if
      end do
    end do

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