function sample!(fin_el::Triangle, s::Matrix{Float64} , wt::Vector{Float64})
  #
  # This subroutine returns the local coordinates and weighting coefficients
  # of the integrating points.
  #
  
  root3 = 1.0 / sqrt(3.0)
  r15 = 0.2 * sqrt(15.0)
  nip = size(s,1)

  w = [5.0/9.0, 8.0/9.0, 5.0/9.0]
  v = [5.0/9.0*w; 8.0/9.0*w; 5.0/9.0*w]
 
  if nip == 1
    s[1,1] = 0.333333333333333
    s[1,2] = 0.333333333333333
    wt[1] = 0.500000000000000
  elseif nip == 3
    s[1,1] = 0.500000000000000
    s[1,2] = 0.500000000000000
    s[2,1] = 0.500000000000000
    s[2,2] = 0.000000000000000
    s[3,1] = 0.000000000000000
    s[3,2] = 0.500000000000000
    wt[1:3] .= 0.333333333333333
    wt[:] .= 0.5*wt
  elseif nip == 4
    s[1,1] = 0.6
    s[1,2] = 0.2
    s[2,1] = 0.2
    s[2,2] = 0.6
    s[3,1] = 0.2
    s[3,2] = 0.2
    s[4,1] = 0.333333333333333
    s[4,2] = 0.333333333333333
    wt[1:3] .= 0.520833333333333
    wt[4] = -0.5625
    wt[:] .= 0.5*wt
  elseif nip == 6
    s[1,1] = 0.816847572980459
    s[1,2] = 0.091576213509771
    s[2,1] = 0.091576213509771
    s[2,2] = 0.816847572980459
    s[3,1] = 0.091576213509771
    s[3,2] = 0.091576213509771
    s[4,1] = 0.108103018168070
    s[4,2] = 0.445948490915965
    s[5,1] = 0.445948490915965
    s[5,2] = 0.108103018168070
    s[6,1] = 0.445948490915965
    s[6,2] = 0.445948490915965
    wt[1:3] .= 0.109951743655322
    wt[4:6] .= 0.223381589678011
    wt[:] .= 0.5*wt
  elseif nip == 7
    s[1,1] = 0.333333333333333
    s[1,2] = 0.333333333333333
    s[2,1] = 0.797426985353087
    s[2,2] = 0.101286507323456
    s[3,1] = 0.101286507323456
    s[3,2] = 0.797426985353087
    s[4,1] = 0.101286507323456
    s[4,2] = 0.101286507323456
    s[5,1] = 0.470142064105115
    s[5,2] = 0.059715871789770
    s[6,1] = 0.059715871789770
    s[6,2] = 0.470142064105115
    s[7,1] = 0.470142064105115
    s[7,2] = 0.470142064105115
    wt[1] = 0.225000000000000
    wt[2:4] .= 0.125939180544827
    wt[5:7] .= 0.132394152788506
    wt[:] .= 0.5*wt
  elseif nip == 12
    s[1,1] = 0.873821971016996
    s[1,2] = 0.063089014491502
    s[2,1] = 0.063089014491502
    s[2,2] = 0.873821971016996
    s[3,1] = 0.063089014491502
    s[3,2] = 0.063089014491502
    s[4,1] = 0.501426509658179
    s[4,2] = 0.249286745170910
    s[5,1] = 0.249286745170910
    s[5,2] = 0.501426509658179
    s[6,1] = 0.249286745170910
    s[6,2] = 0.249286745170910
    s[7,1] = 0.053145049844817
    s[7,2] = 0.310352451033784
    s[8,1] = 0.310352451033784
    s[8,2] = 0.053145049844817
    s[9,1] = 0.053145049844817
    s[9,2] = 0.636502499121398
    s[10,1] = 0.310352451033784
    s[10,2] = 0.636502499121398
    s[11,1] = 0.636502499121398
    s[11,2] = 0.053145049844817
    s[12,1] = 0.636502499121398
    s[12,2] = 0.310352451033784
    wt[1:3] .= 0.050844906370207
    wt[4:6] .= 0.116786275726379
    wt[7:12] .= 0.082851075618374
    wt[:] .= 0.5*wt
  elseif nip == 16
    s[1,1] = 0.333333333333333
    s[1,2] = 0.333333333333333
    s[2,1] = 0.658861384496478
    s[2,2] = 0.170569307751761
    s[3,1] = 0.170569307751761
    s[3,2] = 0.658861384496478
    s[4,1] = 0.170569307751761
    s[4,2] = 0.170569307751761
    s[5,1] = 0.898905543365938
    s[5,2] = 0.050547228317031
    s[6,1] = 0.050547228317031
    s[6,2] = 0.898905543365938
    s[7,1] = 0.050547228317031
    s[7,2] = 0.050547228317031
    s[8,1] = 0.081414823414554
    s[8,2] = 0.459292588292723
    s[9,1] = 0.459292588292723
    s[9,2] = 0.081414823414554
    s[10,1] = 0.459292588292723
    s[10,2] = 0.459292588292723
    s[11,1] = 0.008394777409958
    s[11,2] = 0.263112829634638
    s[12,1] = 0.008394777409958
    s[12,2] = 0.728492392955404
    s[13,1] = 0.263112829634638
    s[13,2] = 0.008394777409958
    s[14,1] = 0.263112829634638
    s[14,2] = 0.728492392955404
    s[15,1] = 0.728492392955404
    s[15,2] = 0.008394777409958
    s[16,1] = 0.728492392955404
    s[16,2] = 0.263112829634638
    wt[1] = 0.144315607677787
    wt[2:4] .= 0.103217370534718
    wt[5:7] .= 0.032458497623198
    wt[8:10] .= 0.095091634267284
    wt[11:16] .= 0.027230314174435
    wt[:] .= 0.5*wt
  else
    println("Wrong number of integrating points for a triangle.")
  end
end
