c***************************************************************
c     subroutine compute_xgJ_triangle(xn,detj)
c     Purpose
c     =======
c     Compute global coordinates of the natural coordinates of
c     the Gauss point for a three-node triangle and the
c     determinant of the Jacobian of the isoparametric mapping
c
c***************************************************************
c
      subroutine compute_xgJ_triangle(xn,detj)
      implicit none
#include "include/params2d.h"
#include "include/gausspt_x.h"
#include "include/gausspt_ref.h" 
c
c     natural coordinates in a triangle: FE shape functions for
c     a three-node triangle (include file "gausspt_ref.h")
c
      real*8 detj 
      real*8 xn(NSD,*)  ! coords. of nodes at vertices of a triangle
      real*8 y12, y23, y31
c
      y12 = xn(2,1) - xn(2,2)
      y23 = xn(2,2) - xn(2,3)
      y31 = xn(2,3) - xn(2,1)
c
c     Gauss point (global coordinates)
c
      xbar = xi1*xn(1,1) + xi2*xn(1,2) + xi3*xn(1,3)
      ybar = xi1*xn(2,1) + xi2*xn(2,2) + xi3*xn(2,3)
c
      detj = (xn(1,1) - xn(1,3))*(xn(2,2) - xn(2,3)) -
     &       (xn(1,2) - xn(1,3))*(xn(2,1) - xn(2,3))
c
      return
      end
c
c***************************************************************
c     subroutine gaussq_triangle(nint,pgauss,wgauss)
c     Purpose
c     =======
c     Gauss quadrature points and weights for triangular element
c     input: nint        - number of quadrature points
c     output: pgauss(3,*) - natural coordinates of Gauss points
c             wgauss(*)   - weights of Gauss points
c
c     Comments
c     ========
c     Rules for nint = 1,3,6,13 are taken from the finite 
c     element book: Cook, Malkus, Plesha, "Concepts and Applns
c     of the Finite Element Method"
c     Rules for nint >= 20 are from the paper: Dunavant, D. A.,
c     "High Degree Efficient Symmetrical Gaussian Quadrature
c     Rules for the Triangle", IJNME, vol. 21, 1129-1148, 1985
c
c***************************************************************
c
      subroutine gaussq_triangle(nint,pgauss,wgauss)
      implicit none
#include "include/params2d.h"
      integer nint
      real*8 pgauss(NODES_PER_TRIANGLE,*), wgauss(*)
      integer i
c
      if (nint .eq. 1) then
c
c     zeroth-order accurate for triangular element
c
         pgauss(1,1) = 1.d0/3.d0
         pgauss(2,1) = pgauss(1,1)
         pgauss(3,1) = pgauss(1,1)
         wgauss(1) = 1.d0
      elseif (nint .eq. 3) then
c
c     second-order accurate for triangular element
c
         pgauss(1,1) = 2.d0/3.d0
         pgauss(2,1) = 1.d0/6.d0
         pgauss(3,1) = pgauss(2,1)
         pgauss(1,2) = pgauss(2,1)
         pgauss(2,2) = pgauss(1,1)
         pgauss(3,2) = pgauss(2,1)
         pgauss(1,3) = pgauss(2,1)
         pgauss(2,3) = pgauss(2,1)
         pgauss(3,3) = pgauss(1,1)
         do i = 1,nint
            wgauss(i) = 1.d0/3.d0
         enddo
c
      elseif (nint .eq. 6) then
c
c     fourth-order accurate for triangular element
c
         pgauss(1,1) = 0.81684 75729 80459d0
         pgauss(2,1) = 0.09157 62135 09771d0
         pgauss(3,1) = pgauss(2,1)
         pgauss(1,2) = pgauss(2,1)
         pgauss(2,2) = pgauss(1,1)
         pgauss(3,2) = pgauss(2,1)
         pgauss(1,3) = pgauss(2,1)
         pgauss(2,3) = pgauss(2,1)
         pgauss(3,3) = pgauss(1,1)
         pgauss(1,4) = 0.10810 30181 68070d0
         pgauss(2,4) = 0.44594 84909 15965d0
         pgauss(3,4) = pgauss(2,4)
         pgauss(1,5) = pgauss(2,4)
         pgauss(2,5) = pgauss(1,4)
         pgauss(3,5) = pgauss(2,4)
         pgauss(1,6) = pgauss(2,4)
         pgauss(2,6) = pgauss(2,4)
         pgauss(3,6) = pgauss(1,4)
         do i = 1,3
            wgauss(i) = 0.10995 17436 55322d0
            wgauss(3+i) = 0.22338 15896 78011d0
         enddo
      elseif (nint .eq. 13) then
c
c     seventh-order accurate for triangular element
c
         pgauss(1,1) = 1.d0/3.d0
         pgauss(2,1) = pgauss(1,1)
         pgauss(3,1) = pgauss(1,1)
c
         pgauss(1,2) = 0.47930 80678 41923d0
         pgauss(2,2) = 0.26034 59660 79038d0
         pgauss(3,2) = pgauss(2,2)
         pgauss(1,3) = pgauss(2,2)
         pgauss(2,3) = pgauss(2,2)
         pgauss(3,3) = pgauss(1,2)
         pgauss(1,4) = pgauss(2,2)
         pgauss(2,4) = pgauss(1,2)
         pgauss(3,4) = pgauss(2,2)
c
         pgauss(1,5) = 0.86973 97941 95568d0
         pgauss(2,5) = 0.06513 01029 02216d0
         pgauss(3,5) = pgauss(2,5)
         pgauss(1,6) = pgauss(2,5)
         pgauss(2,6) = pgauss(2,5)
         pgauss(3,6) = pgauss(1,5)
         pgauss(1,7) = pgauss(2,5)
         pgauss(2,7) = pgauss(1,5)
         pgauss(3,7) = pgauss(2,5)
c
         pgauss(1,8) = 0.63844 41885 69809d0
         pgauss(2,8) = 0.31286 54960 04875d0
         pgauss(3,8) = 0.04869 03154 25316d0
         pgauss(1,9) = pgauss(1,8)
         pgauss(2,9) = pgauss(3,8)
         pgauss(3,9) = pgauss(2,8)
         pgauss(1,10) = pgauss(2,8)
         pgauss(2,10) = pgauss(1,8)
         pgauss(3,10) = pgauss(3,8)
         pgauss(1,11) = pgauss(2,8)
         pgauss(2,11) = pgauss(3,8)
         pgauss(3,11) = pgauss(1,8)
         pgauss(1,12) = pgauss(3,8)
         pgauss(2,12) = pgauss(1,8)
         pgauss(3,12) = pgauss(2,8) 
         pgauss(1,13) = pgauss(3,8)
         pgauss(2,13) = pgauss(2,8)
         pgauss(3,13) = pgauss(1,8)
c
         wgauss(1) = -0.14957 00444 67670d0
         do i = 2,4
            wgauss(i) = 0.17561 52574 33204d0
            wgauss(3+i) = 0.05334 72356 08839d0
         enddo
         do i = 8,13
            wgauss(i) = 0.07711 37608 90257d0
         enddo
c
      elseif (nint .eq. 25) then
c
c     tenth-order accurate for triangular element
c
         pgauss(1,1) = 1.d0/3.d0
         pgauss(2,1) = pgauss(1,1)
         pgauss(3,1) = pgauss(1,1)
c
         pgauss(1,2) = 0.02884 47332 32685d0
         pgauss(2,2) = 0.48557 76333 83657d0
         pgauss(3,2) = pgauss(2,2)
         pgauss(1,3) = pgauss(2,2)
         pgauss(2,3) = pgauss(1,2)
         pgauss(3,3) = pgauss(2,2)
         pgauss(1,4) = pgauss(2,2)
         pgauss(2,4) = pgauss(2,2)
         pgauss(3,4) = pgauss(1,2)
c
         pgauss(1,5) = 0.78103 68490 29926d0
         pgauss(2,5) = 0.10948 15754 85037d0
         pgauss(3,5) = pgauss(2,5)
         pgauss(1,6) = pgauss(2,5)
         pgauss(2,6) = pgauss(1,5)
         pgauss(3,6) = pgauss(2,5)
         pgauss(1,7) = pgauss(2,5)
         pgauss(2,7) = pgauss(2,5)
         pgauss(3,7) = pgauss(1,5)
c
         pgauss(1,8) = 0.14170 72194 14880d0
         pgauss(2,8) = 0.30793 98387 64121d0
         pgauss(3,8) = 0.55035 29418 20999d0
         pgauss(1,9) = pgauss(1,8)
         pgauss(2,9) = pgauss(3,8)
         pgauss(3,9) = pgauss(2,8)
         pgauss(1,10) = pgauss(2,8)
         pgauss(2,10) = pgauss(1,8)
         pgauss(3,10) = pgauss(3,8)
         pgauss(1,11) = pgauss(2,8)
         pgauss(2,11) = pgauss(3,8)
         pgauss(3,11) = pgauss(1,8)
         pgauss(1,12) = pgauss(3,8)
         pgauss(2,12) = pgauss(1,8)
         pgauss(3,12) = pgauss(2,8)
         pgauss(1,13) = pgauss(3,8)
         pgauss(2,13) = pgauss(2,8)
         pgauss(3,13) = pgauss(1,8)
c
         pgauss(1,14) = 0.02500 35347 62686d0
         pgauss(2,14) = 0.24667 25606 39903d0
         pgauss(3,14) = 0.72832 39045 97411d0
         pgauss(1,15) = pgauss(1,14)
         pgauss(2,15) = pgauss(3,14)
         pgauss(3,15) = pgauss(2,14)
         pgauss(1,16) = pgauss(2,14)
         pgauss(2,16) = pgauss(1,14)
         pgauss(3,16) = pgauss(3,14)
         pgauss(1,17) = pgauss(2,14)
         pgauss(2,17) = pgauss(3,14)
         pgauss(3,17) = pgauss(1,14)
         pgauss(1,18) = pgauss(3,14)
         pgauss(2,18) = pgauss(1,14)
         pgauss(3,18) = pgauss(2,14)
         pgauss(1,19) = pgauss(3,14)
         pgauss(2,19) = pgauss(2,14)
         pgauss(3,19) = pgauss(1,14)
c
         pgauss(1,20) = 0.00954 08154 00299d0
         pgauss(2,20) = 0.06680 32510 12200d0
         pgauss(3,20) = 0.92365 59335 87500d0
         pgauss(1,21) = pgauss(1,20)
         pgauss(2,21) = pgauss(3,20)
         pgauss(3,21) = pgauss(2,20)
         pgauss(1,22) = pgauss(2,20)
         pgauss(2,22) = pgauss(1,20)
         pgauss(3,22) = pgauss(3,20)
         pgauss(1,23) = pgauss(2,20)
         pgauss(2,23) = pgauss(3,20)
         pgauss(3,23) = pgauss(1,20)
         pgauss(1,24) = pgauss(3,20)
         pgauss(2,24) = pgauss(1,20)
         pgauss(3,24) = pgauss(2,20)
         pgauss(1,25) = pgauss(3,20)
         pgauss(2,25) = pgauss(2,20)
         pgauss(3,25) = pgauss(1,20)
c
         wgauss(1) = 0.09081 79903 82754d0
         do i = 2,4
            wgauss(i) = 0.03672 59577 56467d0
            wgauss(i+3) = 0.04532 10594 35528d0
         enddo
         do i = 8,13
            wgauss(i)    = 0.07275 79168 45420d0
            wgauss(i+6)  = 0.02832 72425 31057d0
            wgauss(i+12) = 0.00942 16669 63733d0
         enddo
      endif
c
      return
      end
c
c***************************************************************
c     subroutine feshape2d_tri(tri_no,n,ngbrs_list,phi)
c     Purpose
c     =======
c     Compute FE shape functions and derivatives for a 3-node
c     triangular element
c
c***************************************************************
c
      subroutine feshape2d_tri(tri_no,n,ngbrs_list,phi)
      implicit none
#include "include/params2d.h"
#include "include/node.h"
#include "include/triangle.h"
#include "include/gausspt_ref.h" 
c
      integer tri_no, n, ngbrs_list(*)
      real*8 phi(NSD+1,*)
      real*8 xn(NSD,NODES_PER_TRIANGLE)  ! coords. of nodes at
c                                        ! vertices of a triangle
      real*8 x13, x23, y13, y23, detj
c
      integer i
c
      n = 3
      do 10 i = 1,NODES_PER_TRIANGLE
         ngbrs_list(i) = tri_connect(i,tri_no)
         xn(1,i) = coord(1,ngbrs_list(i))
         xn(2,i) = coord(2,ngbrs_list(i))
10    continue  
c
      x13 = xn(1,1) - xn(1,3)
      x23 = xn(1,2) - xn(1,3)
      y13 = xn(2,1) - xn(2,3)
      y23 = xn(2,2) - xn(2,3)     
c
      detj = (xn(1,1) - xn(1,3))*(xn(2,2) - xn(2,3)) -
     &       (xn(1,2) - xn(1,3))*(xn(2,1) - xn(2,3))
c
      phi(1,1) = xi1
      phi(2,1) = y23/detj
      phi(3,1) = -x23/detj
      phi(1,2) = xi2
      phi(2,2) = -y13/detj
      phi(3,2) = x13/detj
      phi(1,3) = xi3
      phi(2,3) = (y13 - y23)/detj
      phi(3,3) = (x23 - x13)/detj
c
      return
      end
