is_in_triangle_first Module Subroutine

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

Arguments

Type IntentOptional Attributes Name
class(type_triangle_first), 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_first~~CallsGraph proc~is_in_triangle_first is_in_triangle_first interface~jacobian_det_triangle_first type_triangle_first%jacobian_det_triangle_first proc~is_in_triangle_first->interface~jacobian_det_triangle_first interface~jacobian_triangle_first type_triangle_first%jacobian_triangle_first proc~is_in_triangle_first->interface~jacobian_triangle_first interface~psi_triangle_first type_triangle_first%psi_triangle_first proc~is_in_triangle_first->interface~psi_triangle_first none~set~3 type_dp_vector_3d%set proc~is_in_triangle_first->none~set~3 proc~get_coordinate abst_mesh%get_coordinate proc~is_in_triangle_first->proc~get_coordinate proc~get_num_nodes~3 abst_mesh%get_num_nodes proc~is_in_triangle_first->proc~get_num_nodes~3 proc~jacobian_det_triangle_first jacobian_det_triangle_first interface~jacobian_det_triangle_first->proc~jacobian_det_triangle_first proc~jacobian_triangle_first jacobian_triangle_first interface~jacobian_triangle_first->proc~jacobian_triangle_first proc~psi_triangle_first psi_triangle_first interface~psi_triangle_first->proc~psi_triangle_first 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_first->interface~jacobian_triangle_first proc~jacobian_triangle_first->proc~get_coordinate proc~jacobian_triangle_first->proc~get_num_nodes~3 interface~dpsi_triangle_first type_triangle_first%dpsi_triangle_first proc~jacobian_triangle_first->interface~dpsi_triangle_first proc~dpsi_triangle_first dpsi_triangle_first interface~dpsi_triangle_first->proc~dpsi_triangle_first

Called by

proc~~is_in_triangle_first~~CalledByGraph proc~is_in_triangle_first is_in_triangle_first interface~is_in_triangle_first type_triangle_first%is_in_triangle_first interface~is_in_triangle_first->proc~is_in_triangle_first

Source Code

    module subroutine is_in_triangle_first(self, cartesian, normalized, is_in)
        implicit none
        class(type_triangle_first), 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_iterations
        real(real64) :: tol
        integer(int32) :: i
        logical :: converged

        call r%set(0.0d0, 0.0d0, 0.0d0)
        tol = 1.0d-15
        max_iterations = 10
        converged = .false.

        ! Newton-Raphson 法による逆写像
        do iter = 1, max_iterations
            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-15) then
                is_in = .false.
                return
            end if

            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 >= -tol) .and. (r%y >= -tol) .and. (r%x + r%y <= 1.0d0 + tol)
        if (is_in) normalized = r

    end subroutine is_in_triangle_first