is_in_triangle_second Module Subroutine

module subroutine is_in_triangle_second(self, cartesian, normalized, is_in)

Arguments

Type IntentOptional Attributes Name
class(type_triangle_second), intent(in) :: self
type(type_dp_vector_3d), intent(in) :: cartesian
type(type_dp_vector_3d), intent(inout) :: normalized
logical, intent(inout) :: is_in

Calls

proc~~is_in_triangle_second~~CallsGraph proc~is_in_triangle_second is_in_triangle_second interface~jacobian_det_triangle_second type_triangle_second%jacobian_det_triangle_second proc~is_in_triangle_second->interface~jacobian_det_triangle_second interface~jacobian_triangle_second type_triangle_second%jacobian_triangle_second proc~is_in_triangle_second->interface~jacobian_triangle_second interface~psi_triangle_second type_triangle_second%psi_triangle_second proc~is_in_triangle_second->interface~psi_triangle_second none~set~3 type_dp_vector_3d%set proc~is_in_triangle_second->none~set~3 proc~get_coordinate abst_mesh%get_coordinate proc~is_in_triangle_second->proc~get_coordinate proc~get_num_nodes~3 abst_mesh%get_num_nodes proc~is_in_triangle_second->proc~get_num_nodes~3 proc~jacobian_det_triangle_second jacobian_det_triangle_second interface~jacobian_det_triangle_second->proc~jacobian_det_triangle_second proc~jacobian_triangle_second jacobian_triangle_second interface~jacobian_triangle_second->proc~jacobian_triangle_second proc~psi_triangle_second psi_triangle_second interface~psi_triangle_second->proc~psi_triangle_second proc~set_dp_vector_3d type_dp_vector_3d%set_dp_vector_3d none~set~3->proc~set_dp_vector_3d proc~set_dp_vector_3d_array type_dp_vector_3d%set_dp_vector_3d_array none~set~3->proc~set_dp_vector_3d_array proc~jacobian_det_triangle_second->interface~jacobian_triangle_second proc~jacobian_triangle_second->proc~get_coordinate proc~jacobian_triangle_second->proc~get_num_nodes~3 interface~dpsi_triangle_second type_triangle_second%dpsi_triangle_second proc~jacobian_triangle_second->interface~dpsi_triangle_second proc~dpsi_triangle_second dpsi_triangle_second interface~dpsi_triangle_second->proc~dpsi_triangle_second

Called by

proc~~is_in_triangle_second~~CalledByGraph proc~is_in_triangle_second is_in_triangle_second interface~is_in_triangle_second type_triangle_second%is_in_triangle_second interface~is_in_triangle_second->proc~is_in_triangle_second

Source Code

    module subroutine is_in_triangle_second(self, cartesian, normalized, is_in)
        class(type_triangle_second), intent(in) :: self
        type(type_dp_vector_3d), intent(in) :: cartesian
        type(type_dp_vector_3d), intent(inout) :: normalized
        logical, intent(inout) :: is_in

        type(type_dp_vector_3d) :: r
        type(type_dp_vector_3d) :: coordinate
        real(real64) :: x0, y0
        real(real64) :: dx_xi, dx_eta, dy_xi, dy_eta
        real(real64) :: detJ
        real(real64) :: dx, dy
        integer(int32) :: iter, max_iter
        real(real64) :: tol
        integer(int32) :: i
        logical :: converged

        ! 初期化
        call r%set(1.0d0 / 3.0d0, 1.0d0 / 3.0d0, 0.0d0)
        tol = 1.0d-10
        max_iter = 30
        converged = .false.

        ! Newton-Raphson 法による逆写像
        do iter = 1, max_iter
            x0 = 0.0d0
            y0 = 0.0d0

            do i = 1, self%get_num_nodes()
                coordinate = self%get_coordinate(i)
                x0 = x0 + self%psi(i, r) * coordinate%x
                y0 = y0 + self%psi(i, r) * coordinate%y
            end do

            dx = cartesian%x - x0
            dy = cartesian%y - y0

            if (sqrt(dx * dx + dy * dy) < tol) then
                converged = .true.
                exit
            end if

            dx_xi = self%jacobian(1, 1, r)
            dx_eta = self%jacobian(1, 2, r)
            dy_xi = self%jacobian(2, 1, r)
            dy_eta = self%jacobian(2, 2, r)

            detJ = self%jacobian_Det(r)
            if (abs(detJ) < 1.0d-20) exit ! ヤコビ行列の特異性チェック

            ! Newton-Raphson 更新
            r%x = r%x + (dy_eta * dx - dx_eta * dy) / detJ
            r%y = r%y + (-dy_xi * dx + dx_xi * dy) / detJ
        end do

        ! 最終判定：収束かつ自然座標が範囲内
        is_in = converged .and. (r%x >= 0.0d0) .and. (r%y >= 0.0d0) .and. (r%x + r%y <= 1.0d0)
        if (is_in) normalized = r

    end subroutine is_in_triangle_second