Commits

Lisandro Dalcin committed 707370f

Add Fortran implementation of quadrature rules

Comments (0)

Files changed (2)

 SOURCEH  = ../include/petiga.h petigapc.h petigagrid.h petigapart.h
 SOURCEC  = petiga.c petigareg.c petigaaxis.c petigarule.c petigabasis.c petigabound.c petigaelem.c petigapoint.c petigavec.c petigamat.c petigascl.c petigapcb.c petigapce.c petigaksp.c petigasnes.c petigats.c petigats2.c petigaio.c petigagrid.c petigapart.c snesfdcolor.c tsalpha2.c
 SOURCEF1 = petigaftn.F90 petigaval.F90
-SOURCEF2 = petigabsp.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90
+SOURCEF2 = petigabsp.f90 petigaqdr.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90
 SOURCEF  = ${SOURCEF1} ${SOURCEF2}
 OBJSC    = ${SOURCEC:.c=.o}
 OBJSF    = ${SOURCEF1:.F90=.o} ${SOURCEF2:.f90=.o}

src/petigaqdr.f90

+! -*- 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
+     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