Source

PetIGA / src / petigaqdr.f90

Full commit
! -*- f90 -*-

subroutine IGA_Rule_GaussLegendre(q,X,W) &
  bind(C, name="IGA_Rule_GaussLegendre")
  use PetIGA
  implicit none
  integer(kind=IGA_INTEGER_KIND), intent(in)  :: q
  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(0:q-1)
  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(0:q-1)
  select case (q)
  case (1) ! p <= 1
     X(0) = 0.0
     W(0) = 2.0
  case (2) ! p <= 3
     X(0) = -0.5773502691896257645091487805019576 ! 1/sqrt(3)
     X(1) = -X(0)
     W(0) =  1.0
     W(1) =  W(0)
  case (3) ! p <= 5
     X(0) = -0.7745966692414833770358530799564799 ! sqrt(3/5)
     X(1) =  0.0
     X(2) = -X(0)
     W(0) =  0.5555555555555555555555555555555556 ! 5/9
     W(1) =  0.8888888888888888888888888888888889 ! 8/9
     W(2) =  W(0)
  case (4) ! p <= 7
     X(0) = -0.8611363115940525752239464888928094 ! sqrt((3+2*sqrt(6/5))/7)
     X(1) = -0.3399810435848562648026657591032448 ! sqrt((3-2*sqrt(6/5))/7)
     X(2) = -X(1)
     X(3) = -X(0)
     W(0) =  0.3478548451374538573730639492219994 ! (18-sqrt(30))/36
     W(1) =  0.6521451548625461426269360507780006 ! (18+sqrt(30))/36
     W(2) =  W(1)
     W(3) =  W(0)
  case (5) ! p <= 9
     X(0) = -0.9061798459386639927976268782993929 ! 1/3*sqrt(5+2*sqrt(10/7))
     X(1) = -0.5384693101056830910363144207002086 ! 1/3*sqrt(5-2*sqrt(10/7))
     X(2) =  0.0
     X(3) = -X(1)
     X(4) = -X(0)
     W(0) =  0.2369268850561890875142640407199173 ! (322-13*sqrt(70))/900
     W(1) =  0.4786286704993664680412915148356382 ! (322+13*sqrt(70))/900
     W(2) =  0.5688888888888888888888888888888889 ! 128/225
     W(3) =  W(1)
     W(4) =  W(0)
  case (6) ! p <= 11
     X(0) = -0.9324695142031520278123015544939946 ! << NumericalDifferentialEquationAnalysis`
     X(1) = -0.6612093864662645136613995950199053 ! GaussianQuadratureWeights(6, -1, 1, 37)
     X(2) = -0.2386191860831969086305017216807119
     X(3) = -X(2)
     X(4) = -X(1)
     X(5) = -X(0)
     W(0) =  0.171324492379170345040296142172732894
     W(1) =  0.360761573048138607569833513837716112
     W(2) =  0.467913934572691047389870343989550995
     W(3) =  W(2)
     W(4) =  W(1)
     W(5) =  W(0)
  case (7) ! p <= 13
     X(0) = -0.9491079123427585245261896840478513 ! << NumericalDifferentialEquationAnalysis`
     X(1) = -0.7415311855993944398638647732807884 ! GaussianQuadratureWeights(7, -1, 1, 37)
     X(2) = -0.4058451513773971669066064120769615
     X(3) =  0.0
     X(4) = -X(2)
     X(5) = -X(1)
     X(6) = -X(0)
     W(0) =  0.129484966168869693270611432679082018
     W(1) =  0.279705391489276667901467771423779582
     W(2) =  0.381830050505118944950369775488975134
     W(3) =  0.417959183673469387755102040816326531
     W(4) =  W(2)
     W(5) =  W(1)
     W(6) =  W(0)
  case (8) ! p <= 15
     X(0) = -0.9602898564975362316835608685694730 ! << NumericalDifferentialEquationAnalysis`
     X(1) = -0.7966664774136267395915539364758304 ! GaussianQuadratureWeights(8, -1, 1, 37)
     X(2) = -0.5255324099163289858177390491892463
     X(3) = -0.1834346424956498049394761423601840
     X(4) = -X(3)
     X(5) = -X(2)
     X(6) = -X(1)
     X(7) = -X(0)
     W(0) =  0.101228536290376259152531354309962190
     W(1) =  0.222381034453374470544355994426240884
     W(2) =  0.313706645877887287337962201986601313
     W(3) =  0.362683783378361982965150449277195612
     W(4) =  W(3)
     W(5) =  W(2)
     W(6) =  W(1)
     W(7) =  W(0)
  case (9) ! p <= 17
     X(0) = -0.9681602395076260898355762029036729 ! << NumericalDifferentialEquationAnalysis`
     X(1) = -0.8360311073266357942994297880697349 ! GaussianQuadratureWeights(9, -1, 1, 37)
     X(2) = -0.6133714327005903973087020393414742
     X(3) = -0.3242534234038089290385380146433366
     X(4) =  0.0
     X(5) = -X(3)
     X(6) = -X(2)
     X(7) = -X(1)
     X(8) = -X(0)
     W(0) =  0.081274388361574411971892158110523651
     W(1) =  0.180648160694857404058472031242912810
     W(2) =  0.260610696402935462318742869418632850
     W(3) =  0.312347077040002840068630406584443666
     W(4) =  0.330239355001259763164525069286974049
     W(5) =  W(3)
     W(6) =  W(2)
     W(7) =  W(1)
     W(8) =  W(0)
  case default
     X = 0.0
     W = 0.0
  end select
end subroutine IGA_Rule_GaussLegendre

subroutine IGA_Rule_GaussLobatto(q,X,W) &
  bind(C, name="IGA_Rule_GaussLobatto")
  use PetIGA
  implicit none
  integer(kind=IGA_INTEGER_KIND), intent(in)  :: q
  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(0:q-1)
  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(0:q-1)
  select case (q)
  case (2) ! p <= 1
     X(0) = -1.0
     X(1) = -X(0)
     W(0) =  1.0
     W(1) =  W(0)
  case (3) ! p <= 3
     X(0) = -1.0
     X(1) =  0.0
     X(2) = -X(0)
     W(0) =  0.3333333333333333333333333333333333  ! 1/3
     W(1) =  1.3333333333333333333333333333333333  ! 4/3
     W(2) =  W(0)
  case (4) ! p <= 5
     X(0) = -1.0
     X(1) = -0.44721359549995793928183473374625525 ! 1/sqrt(5)
     X(2) = -X(1)
     X(3) = -X(0)
     W(0) =  0.16666666666666666666666666666666667 ! 1/6
     W(1) =  0.83333333333333333333333333333333333 ! 5/6
     W(2) =  W(1)
     W(3) =  W(0)
  case (5) ! p <= 7
     X(0) = -1.0
     X(1) = -0.65465367070797714379829245624685835 ! sqrt(3/7)
     X(2) =  0.0
     X(3) = -X(1)
     X(4) = -X(0)
     W(0) =  0.10000000000000000000000000000000000 !  1/10
     W(1) =  0.54444444444444444444444444444444444 ! 49/90
     W(2) =  0.71111111111111111111111111111111111 ! 32/45
     W(3) =  W(1)
     W(4) =  W(0)
  case default
     X = 0.0
     W = 0.0
  end select
end subroutine IGA_Rule_GaussLobatto