!************************************************************************ !* Educationally-Designed Unstructured 2D (EDU2D) Code !* !* --- EDU2D_Quadrature_Triangle !* !* This program illustrates the use of a Gaussian quadrature, which !* integrates a function over a triangle. !* !* !* Note: Integral over a triangle may be used for computing !* a cell-averaged solution (or errors) on unstructured grids. !* !* !* This is Version 1 (February 22, 2017). !* !* This F90 code is written and made available for an educational purpose. !* This file may be updated in future. !* !* Katate Masatsuka, http://www.cfdbooks.com !************************************************************************ program edu2d_quadrature_triangle implicit none integer , parameter :: p2 = selected_real_kind(p=15) real(p2) :: x1, x2, x3, y1, y2, y3 real(p2) :: integral x1 = 0.1_p2 x2 = 1.3_p2 x3 = -0.9_p2 y1 = 0.3_p2 y2 = -0.2_p2 y3 = 1.7_p2 ! Compute the integral of 'func' by a seven-point quadrature: call volume_quintic(x1,x2,x3, y1,y2,y3, func, integral) write(*,*) " ---------------------------------------------------------- " write(*,*) write(*,'(a)') " Vertex coordinates of a triangle T: " write(*,*) write(*,*) " x1 = ", x1, " y1 = ", y1 write(*,*) " x2 = ", x2, " y2 = ", y2 write(*,*) " x3 = ", x3, " y3 = ", y3 write(*,*) write(*,'(a)') " Integral of 'func' over T: " write(*,*) write(*,'(a)') " Example function = x^5 + y^5 -32*x^2*y^3 + x^3*y^2" write(*,*) write(*,*) write(*,'(a,f25.15)') " 7-point quadrature = ", integral write(*,'(a,f25.15)') " Exact integral = ", 0.53726127666666666668_p2 write(*,*) write(*,*) " (Note: 7-point quadrature is exact for polynomials of degree 5.)" write(*,*) contains !************************************************************************ ! Function to be integrated. ! ! Exmaple: x^5 + y^5 -32*x^2*y^3 + x^3*y^2 ! ! Integral of this polynomial over a triangle defined by ! the vertices: ! (0.1,0.3), (1.3,-0.2), (-0.9,1.7) ! ! is 0.53726127666666666668. ! !************************************************************************ pure function func(x,y) implicit none integer , parameter :: p2 = selected_real_kind(p=15) real(p2), intent(in) :: x, y real(p2) :: func func = x**5 + y**5 + x**2*y**3 + x**3*y**2 end function func !************************************************************************ !************************************************************************ ! ! Seven-point quadrature formula. ! !------------------------------------------------------------------------ ! Input: ! ! (x1,y1), (x2,y2), (x3,y3): Vertex coordinates of a triangle ! on which a function is integrated. ! f(x,y): Function to be integrated. ! Output: ! ! integral : Integral value computed by a 7-point quadrature formula. !------------------------------------------------------------------------ ! ! - This is exact for polynomials of degree 5 (quintic). ! - Error = O(h^6) ! ! Ref.["High Degree Efficint Symmetrical Gaussian Quadrature Rules for ! the Triangle" by D. A. Dunavant, International Journal for ! Numerical Methods in Engineering, Vol.21, 1129-1148, 1985] ! !* Katate Masatsuka, http://www.cfdbooks.com !************************************************************************ subroutine volume_quintic(x1,x2,x3, y1,y2,y3, f, integral) implicit none integer , parameter :: p2 = selected_real_kind(p=15) !Input !(1)Function to be integrated. interface function f(x,y) integer , parameter :: p2 = selected_real_kind(p=15) real(p2), intent(in) :: x, y real(p2) :: f end function end interface !(2)Vertex corrdinates of the triangle on which functn is integrated: real(p2), intent( in) :: x1,x2,x3, y1,y2,y3 !Output real(p2), intent(out) :: integral !Local variables real(p2) :: a, b, c, d, vol, w1, w2, w3 real(p2), dimension(7) :: weight, xc, yc integer :: i ! Volume (area) of the triangle vol = 0.5_p2*( x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2) ) ! Control points for the seven-point Gaussian quadrature a = 0.059715871789770_p2 b = 0.470142064105115_p2 c = 0.797426985353087_p2 d = 0.101286507323456_p2 w1 = 0.225000000000000_p2 w2 = 0.132394152788506_p2 w3 = 0.125939180544827_p2 xc(1) = (x1+x2+x3)/3.0_p2 xc(2) = a*x1 + b*(x2+x3) xc(3) = a*x2 + b*(x1+x3) xc(4) = b*(x1+x2) + a*x3 xc(5) = c*x1 + d*(x2+x3) xc(6) = c*x2 + d*(x1+x3) xc(7) = d*(x1+x2) + c*x3 yc(1) = (y1+y2+y3)/3.0_p2 yc(2) = a*y1 + b*(y2+y3) yc(3) = a*y2 + b*(y1+y3) yc(4) = b*(y1+y2) + a*y3 yc(5) = c*y1 + d*(y2+y3) yc(6) = c*y2 + d*(y1+y3) yc(7) = d*(y1+y2) + c*y3 ! Weights for the seven-point Gaussian quadrature. weight(1) = w1 weight(2) = w2 weight(3) = w2 weight(4) = w2 weight(5) = w3 weight(6) = w3 weight(7) = w3 ! Add up the contributions to integrate. integral = 0.0_p2 loop_control_pts : do i = 1, 7 integral = integral + weight(i) * f( xc(i), yc(i) ) end do loop_control_pts ! Multiply the result by the volume. integral = integral * vol return end subroutine volume_quintic !************************************************************************ end program edu2d_quadrature_triangle