Commits

Lisandro Dalcin committed 886135d

Rename Fortran sources to uppercase extension

Comments (0)

Files changed (12)

 set (CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
 
 file (GLOB PetIGA_SOURCES_C RELATIVE ${PetIGA_SOURCE_DIR} ${PetIGA_SOURCE_DIR}/src/*.c)
-file (GLOB PetIGA_SOURCES_F RELATIVE ${PetIGA_SOURCE_DIR} ${PetIGA_SOURCE_DIR}/src/*.[Ff]90)
+file (GLOB PetIGA_SOURCES_F RELATIVE ${PetIGA_SOURCE_DIR} ${PetIGA_SOURCE_DIR}/src/*.F90)
 set  (PetIGA_SOURCES_ALL ${PetIGA_SOURCES_C} ${PetIGA_SOURCES_F})
 if (PETSC_CLANGUAGE_Cxx)
   foreach (file ${PetIGA_SOURCES_C})
 
 SOURCEH  = ../include/petiga.h petigabl.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 petigapc.c petigaksp.c petigasnes.c petigats.c petigadm.c petigats2.c petigaio.c petigagrid.c petigapart.c petscvwopt.c snesfdcolor.c tsalpha2.c
-SOURCEF1 = petigaftn.F90 petigaval.F90
-SOURCEF2 = petigabsp.f90 petigaqdr.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90
-SOURCEF  = ${SOURCEF1} ${SOURCEF2}
+SOURCEF = petigaftn.F90 petigaval.F90 petigabsp.F90 petigaqdr.F90 petiga1d.F90 petiga2d.F90 petiga3d.F90
 OBJSC    = ${SOURCEC:.c=.o}
-OBJSF    = ${SOURCEF1:.F90=.o} ${SOURCEF2:.f90=.o}
+OBJSF    = ${SOURCEF:.F90=.o}
 
 LIBBASE  = libpetiga
 DIRS     =
 include ${PETIGA_DIR}/conf/petigavariables
 include ${PETIGA_DIR}/conf/petigarules
 include ${PETIGA_DIR}/conf/petigatest
-
-OBJSC    = ${SOURCEC:.c=.o}
-OBJSF    = ${SOURCEF1:.F90=.o} ${SOURCEF2:.f90=.o}
+pure subroutine IGA_Quadrature_1D(&
+     inq,iX,iW,iJ,                &
+     X,W,J)                       &
+  bind(C, name="IGA_Quadrature_1D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iX(inq), iW(inq), iJ
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(    inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(    inq)
+  integer(kind=IGA_INTEGER_KIND)  :: iq
+  do iq=1,inq
+     X(1,iq) = iX(iq)
+     W(  iq) = iW(iq)
+  end do
+  J = iJ
+end subroutine IGA_Quadrature_1D
+
+
+pure subroutine IGA_BasisFuns_1D(&
+     order,                      &
+     rational,W,                 &
+     inq,ina,ind,iN,             &
+     N0,N1,N2,N3)                &
+  bind(C, name="IGA_BasisFuns_1D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: rational
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq, ina, ind
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: W(ina)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(       ina,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(   dim,ina,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(dim**2,ina,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim**3,ina,inq)
+  integer(kind=IGA_INTEGER_KIND)  :: ia, iq
+  integer(kind=IGA_INTEGER_KIND)  :: nen
+  nen = ina
+  do iq=1,inq
+     call TensorBasisFuns(&
+          order,&
+          ina,ind,iN(:,:,iq),&
+          N0(  :,iq),&
+          N1(:,:,iq),&
+          N2(:,:,iq),&
+          N3(:,:,iq))
+     if (rational /= 0) then
+        call Rationalize(&
+             order,&
+             nen,W,&
+             N0(  :,iq),&
+             N1(:,:,iq),&
+             N2(:,:,iq),&
+             N3(:,:,iq))
+     end if
+  end do
+
+contains
+
+pure subroutine TensorBasisFuns(&
+     ord,&
+     ina,ind,iN,&
+     N0,N1,N2,N3)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ord
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ina, ind
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(            ina)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(        dim,ina)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(    dim,dim,ina)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim,dim,dim,ina)
+  integer(kind=IGA_INTEGER_KIND)  :: ia
+  !
+  do ia=1,ina
+     N0(ia) = iN(0,ia)
+  end do
+  !
+  do ia=1,ina
+     N1(1,ia) = iN(1,ia)
+  end do
+  !
+  if (ord < 2) return
+  do ia=1,ina
+     N2(1,1,ia) = iN(2,ia)
+  end do
+  !
+  if (ord < 3) return
+  do ia=1,ina
+     N3(1,1,1,ia) = iN(3,ia)
+  end do
+  !
+end subroutine TensorBasisFuns
+
+include 'petigarat.f90.in'
+
+end subroutine IGA_BasisFuns_1D
+
+
+pure subroutine IGA_ShapeFuns_1D(&
+     order,                      &
+     nqp,nen,X,                  &
+     M0,M1,M2,M3,                &
+     N0,N1,N2,N3,                &
+     dX,G0,G1,H0,H1)             &
+  bind(C, name="IGA_ShapeFuns_1D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nqp
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: X(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M0(       nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M1(dim,   nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M2(dim**2,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M3(dim**3,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N0(       nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
+  integer(kind=IGA_INTEGER_KIND)  :: q
+  do q=1,nqp
+     call GeometryMap(&
+          order,&
+          nen,X,&
+          M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
+          N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
+          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
+  end do
+contains
+include 'petigageo.f90.in'
+end subroutine IGA_ShapeFuns_1D

src/petiga1d.f90

-pure subroutine IGA_Quadrature_1D(&
-     inq,iX,iW,iJ,                &
-     X,W,J)                       &
-  bind(C, name="IGA_Quadrature_1D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iX(inq), iW(inq), iJ
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(    inq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(    inq)
-  integer(kind=IGA_INTEGER_KIND)  :: iq
-  do iq=1,inq
-     X(1,iq) = iX(iq)
-     W(  iq) = iW(iq)
-  end do
-  J = iJ
-end subroutine IGA_Quadrature_1D
-
-
-pure subroutine IGA_BasisFuns_1D(&
-     order,                      &
-     rational,W,                 &
-     inq,ina,ind,iN,             &
-     N0,N1,N2,N3)                &
-  bind(C, name="IGA_BasisFuns_1D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: rational
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq, ina, ind
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: W(ina)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina,inq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(       ina,inq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(   dim,ina,inq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(dim**2,ina,inq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim**3,ina,inq)
-  integer(kind=IGA_INTEGER_KIND)  :: ia, iq
-  integer(kind=IGA_INTEGER_KIND)  :: nen
-  nen = ina
-  do iq=1,inq
-     call TensorBasisFuns(&
-          order,&
-          ina,ind,iN(:,:,iq),&
-          N0(  :,iq),&
-          N1(:,:,iq),&
-          N2(:,:,iq),&
-          N3(:,:,iq))
-     if (rational /= 0) then
-        call Rationalize(&
-             order,&
-             nen,W,&
-             N0(  :,iq),&
-             N1(:,:,iq),&
-             N2(:,:,iq),&
-             N3(:,:,iq))
-     end if
-  end do
-
-contains
-
-pure subroutine TensorBasisFuns(&
-     ord,&
-     ina,ind,iN,&
-     N0,N1,N2,N3)
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ord
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ina, ind
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(            ina)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(        dim,ina)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(    dim,dim,ina)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim,dim,dim,ina)
-  integer(kind=IGA_INTEGER_KIND)  :: ia
-  !
-  do ia=1,ina
-     N0(ia) = iN(0,ia)
-  end do
-  !
-  do ia=1,ina
-     N1(1,ia) = iN(1,ia)
-  end do
-  !
-  if (ord < 2) return
-  do ia=1,ina
-     N2(1,1,ia) = iN(2,ia)
-  end do
-  !
-  if (ord < 3) return
-  do ia=1,ina
-     N3(1,1,1,ia) = iN(3,ia)
-  end do
-  !
-end subroutine TensorBasisFuns
-
-include 'petigarat.f90.in'
-
-end subroutine IGA_BasisFuns_1D
-
-
-pure subroutine IGA_ShapeFuns_1D(&
-     order,                      &
-     nqp,nen,X,                  &
-     M0,M1,M2,M3,                &
-     N0,N1,N2,N3,                &
-     dX,G0,G1,H0,H1)             &
-  bind(C, name="IGA_ShapeFuns_1D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nqp
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: X(dim,nen)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M0(       nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M1(dim,   nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N0(       nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
-  integer(kind=IGA_INTEGER_KIND)  :: q
-  do q=1,nqp
-     call GeometryMap(&
-          order,&
-          nen,X,&
-          M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
-          N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
-  end do
-contains
-include 'petigageo.f90.in'
-end subroutine IGA_ShapeFuns_1D
+pure subroutine IGA_Quadrature_2D(&
+     inq,iX,iW,iJ,                &
+     jnq,jX,jW,jJ,                &
+     X,W,J)                       &
+  bind(C, name="IGA_Quadrature_2D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jnq
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iX(inq), iW(inq), iJ
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jX(jnq), jW(jnq), jJ
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(    inq,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(    inq,jnq)
+  integer(kind=IGA_INTEGER_KIND)  :: iq
+  integer(kind=IGA_INTEGER_KIND)  :: jq
+  do jq=1,jnq; do iq=1,inq
+     X(1,iq,jq) = iX(iq)
+     X(2,iq,jq) = jX(jq)
+     W(  iq,jq) = iW(iq) * jW(jq)
+  end do; end do
+  J = iJ * jJ
+end subroutine IGA_Quadrature_2D
+
+
+pure subroutine IGA_BasisFuns_2D(&
+     order,                      &
+     rational,W,                 &
+     inq,ina,ind,iN,             &
+     jnq,jna,jnd,jN,             &
+     N0,N1,N2,N3)                &
+  bind(C, name="IGA_BasisFuns_2D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: rational
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq, ina, ind
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jnq, jna, jnd
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: W(ina*jna)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(       ina*jna,inq,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(   dim,ina*jna,inq,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(dim**2,ina*jna,inq,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim**3,ina*jna,inq,jnq)
+  integer(kind=IGA_INTEGER_KIND)  :: ia, iq
+  integer(kind=IGA_INTEGER_KIND)  :: ja, jq
+  integer(kind=IGA_INTEGER_KIND)  :: ka, kq
+  integer(kind=IGA_INTEGER_KIND)  :: nen
+  nen = ina*jna
+  do jq=1,jnq
+     do iq=1,inq
+        call TensorBasisFuns(&
+             order,&
+             ina,ind,iN(:,:,iq),&
+             jna,jnd,jN(:,:,jq),&
+             N0(  :,iq,jq),&
+             N1(:,:,iq,jq),&
+             N2(:,:,iq,jq),&
+             N3(:,:,iq,jq))
+        if (rational /= 0) then
+           call Rationalize(&
+                order,&
+                nen,W,&
+                N0(  :,iq,jq),&
+                N1(:,:,iq,jq),&
+                N2(:,:,iq,jq),&
+                N3(:,:,iq,jq))
+        end if
+     end do
+  end do
+
+contains
+
+pure subroutine TensorBasisFuns(&
+     ord,&
+     ina,ind,iN,&
+     jna,jnd,jN,&
+     N0,N1,N2,N3)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ord
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ina, ind
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jna, jnd
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(            ina,jna)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(        dim,ina,jna)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(    dim,dim,ina,jna)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim,dim,dim,ina,jna)
+  integer(kind=IGA_INTEGER_KIND)  :: ia, ja
+  !
+   do ja=1,jna; do ia=1,ina
+     N0(ia,ja) = iN(0,ia) * jN(0,ja)
+  end do; end do
+  !
+  do ja=1,jna; do ia=1,ina
+     N1(1,ia,ja) = iN(1,ia) * jN(0,ja)
+     N1(2,ia,ja) = iN(0,ia) * jN(1,ja)
+  end do; end do
+  !
+  if (ord < 2) return ! XXX Optimize!
+  do ja=1,jna; do ia=1,ina
+     N2(1,1,ia,ja) = iN(2,ia) * jN(0,ja)
+     N2(2,1,ia,ja) = iN(1,ia) * jN(1,ja)
+     N2(1,2,ia,ja) = iN(1,ia) * jN(1,ja)
+     N2(2,2,ia,ja) = iN(0,ia) * jN(2,ja)
+  end do; end do
+  !
+  if (ord < 3) return ! XXX Optimize!
+  do ja=1,jna; do ia=1,ina
+     N3(1,1,1,ia,ja) = iN(3,ia) * jN(0,ja)
+     N3(2,1,1,ia,ja) = iN(2,ia) * jN(1,ja)
+     N3(1,2,1,ia,ja) = iN(2,ia) * jN(1,ja)
+     N3(2,2,1,ia,ja) = iN(1,ia) * jN(2,ja)
+     N3(1,1,2,ia,ja) = iN(2,ia) * jN(1,ja)
+     N3(2,1,2,ia,ja) = iN(1,ia) * jN(2,ja)
+     N3(1,2,2,ia,ja) = iN(1,ia) * jN(2,ja)
+     N3(2,2,2,ia,ja) = iN(0,ia) * jN(3,ja)
+  end do; end do
+  !
+end subroutine TensorBasisFuns
+
+include 'petigarat.f90.in'
+
+end subroutine IGA_BasisFuns_2D
+
+
+pure subroutine IGA_ShapeFuns_2D(&
+     order,                      &
+     nqp,nen,X,                  &
+     M0,M1,M2,M3,                &
+     N0,N1,N2,N3,                &
+     dX,G0,G1,H0,H1)             &
+  bind(C, name="IGA_ShapeFuns_2D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nqp
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: X(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M0(       nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M1(dim,   nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M2(dim**2,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M3(dim**3,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N0(       nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
+  integer(kind=IGA_INTEGER_KIND)  :: q
+  do q=1,nqp
+     call GeometryMap(&
+          order,&
+          nen,X,&
+          M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
+          N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
+          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
+  end do
+contains
+include 'petigageo.f90.in'
+end subroutine IGA_ShapeFuns_2D
+
+
+subroutine IGA_BoundaryArea_2D(&
+     m,axis,side,              &
+     geometry,Cx,              &
+     rational,Cw,              &
+     nqp,W,nen,ndr,N,          &
+     dS)                       &
+  bind(C, name="IGA_BoundaryArea_2D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter         :: nsd = 2
+  integer(kind=IGA_INTEGER_KIND), parameter         :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in)        :: m(2)
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: axis, side
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: geometry, rational
+  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cx(nsd,m(1),m(2))
+  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cw(    m(1),m(2))
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: nqp, nen, ndr
+  real   (kind=IGA_REAL_KIND   ), intent(in)        :: W(nqp), N(0:ndr,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)       :: dS
+  integer(kind=IGA_INTEGER_KIND)  :: k, q
+  real   (kind=IGA_REAL_KIND   )  :: N0(nen), N1(dim,nen), detJ
+  real   (kind=IGA_REAL_KIND   ), pointer :: Xx(:,:), Xw(:)
+  select case (axis)
+  case (0)
+     if (side==0) k=1
+     if (side==1) k=m(1)
+     Xx => Cx(:,k,:); Xw => Cw(k,:)
+  case (1)
+     if (side==0) k=1
+     if (side==1) k=m(2)
+     Xx => Cx(:,:,k); Xw => Cw(:,k)
+  end select
+  detJ = 1
+  dS = 0
+  do q=1,nqp
+     N0(  :) = N(0,:,q)
+     N1(1,:) = N(1,:,q)
+     if (rational /= 0) then
+        call Rationalize(nen,Xw,N0,N1)
+     end if
+     if (geometry /= 0) then
+        call Jacobian(nen,N1,Xx,detJ)
+     end if
+     dS = dS + detJ * W(q)
+  end do
+contains
+pure subroutine Rationalize(nen,W,R0,R1)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter     :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in)    :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: W(nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R0(    nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R1(dim,nen)
+  integer(kind=IGA_INTEGER_KIND)  :: i
+  real   (kind=IGA_REAL_KIND   )  :: W0
+  R0 = W * R0
+  W0 = sum(R0)
+  R0 = R0 / W0
+  do i=1,dim
+  R1(i,:) = W*R1(i,:) - R0 * sum(W*R1(i,:))
+  end do
+  R1 = R1 / W0
+end subroutine Rationalize
+pure subroutine Jacobian(nen,N,X,J)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: nsd = 2
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: X(nsd,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: J
+  real   (kind=IGA_REAL_KIND   )  :: F(dim,nsd), M(dim,dim)
+  F = matmul(N,transpose(X))
+  M = matmul(F,transpose(F))
+  J = sqrt(abs(Determinant(dim,M)))
+end subroutine Jacobian
+include 'petigainv.f90.in'
+end subroutine IGA_BoundaryArea_2D

src/petiga2d.f90

-pure subroutine IGA_Quadrature_2D(&
-     inq,iX,iW,iJ,                &
-     jnq,jX,jW,jJ,                &
-     X,W,J)                       &
-  bind(C, name="IGA_Quadrature_2D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jnq
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iX(inq), iW(inq), iJ
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jX(jnq), jW(jnq), jJ
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq,jnq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(    inq,jnq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(    inq,jnq)
-  integer(kind=IGA_INTEGER_KIND)  :: iq
-  integer(kind=IGA_INTEGER_KIND)  :: jq
-  do jq=1,jnq; do iq=1,inq
-     X(1,iq,jq) = iX(iq)
-     X(2,iq,jq) = jX(jq)
-     W(  iq,jq) = iW(iq) * jW(jq)
-  end do; end do
-  J = iJ * jJ
-end subroutine IGA_Quadrature_2D
-
-
-pure subroutine IGA_BasisFuns_2D(&
-     order,                      &
-     rational,W,                 &
-     inq,ina,ind,iN,             &
-     jnq,jna,jnd,jN,             &
-     N0,N1,N2,N3)                &
-  bind(C, name="IGA_BasisFuns_2D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: rational
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq, ina, ind
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jnq, jna, jnd
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: W(ina*jna)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina,inq)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna,jnq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(       ina*jna,inq,jnq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(   dim,ina*jna,inq,jnq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(dim**2,ina*jna,inq,jnq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim**3,ina*jna,inq,jnq)
-  integer(kind=IGA_INTEGER_KIND)  :: ia, iq
-  integer(kind=IGA_INTEGER_KIND)  :: ja, jq
-  integer(kind=IGA_INTEGER_KIND)  :: ka, kq
-  integer(kind=IGA_INTEGER_KIND)  :: nen
-  nen = ina*jna
-  do jq=1,jnq
-     do iq=1,inq
-        call TensorBasisFuns(&
-             order,&
-             ina,ind,iN(:,:,iq),&
-             jna,jnd,jN(:,:,jq),&
-             N0(  :,iq,jq),&
-             N1(:,:,iq,jq),&
-             N2(:,:,iq,jq),&
-             N3(:,:,iq,jq))
-        if (rational /= 0) then
-           call Rationalize(&
-                order,&
-                nen,W,&
-                N0(  :,iq,jq),&
-                N1(:,:,iq,jq),&
-                N2(:,:,iq,jq),&
-                N3(:,:,iq,jq))
-        end if
-     end do
-  end do
-
-contains
-
-pure subroutine TensorBasisFuns(&
-     ord,&
-     ina,ind,iN,&
-     jna,jnd,jN,&
-     N0,N1,N2,N3)
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ord
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ina, ind
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jna, jnd
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(            ina,jna)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(        dim,ina,jna)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(    dim,dim,ina,jna)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim,dim,dim,ina,jna)
-  integer(kind=IGA_INTEGER_KIND)  :: ia, ja
-  !
-   do ja=1,jna; do ia=1,ina
-     N0(ia,ja) = iN(0,ia) * jN(0,ja)
-  end do; end do
-  !
-  do ja=1,jna; do ia=1,ina
-     N1(1,ia,ja) = iN(1,ia) * jN(0,ja)
-     N1(2,ia,ja) = iN(0,ia) * jN(1,ja)
-  end do; end do
-  !
-  if (ord < 2) return ! XXX Optimize!
-  do ja=1,jna; do ia=1,ina
-     N2(1,1,ia,ja) = iN(2,ia) * jN(0,ja)
-     N2(2,1,ia,ja) = iN(1,ia) * jN(1,ja)
-     N2(1,2,ia,ja) = iN(1,ia) * jN(1,ja)
-     N2(2,2,ia,ja) = iN(0,ia) * jN(2,ja)
-  end do; end do
-  !
-  if (ord < 3) return ! XXX Optimize!
-  do ja=1,jna; do ia=1,ina
-     N3(1,1,1,ia,ja) = iN(3,ia) * jN(0,ja)
-     N3(2,1,1,ia,ja) = iN(2,ia) * jN(1,ja)
-     N3(1,2,1,ia,ja) = iN(2,ia) * jN(1,ja)
-     N3(2,2,1,ia,ja) = iN(1,ia) * jN(2,ja)
-     N3(1,1,2,ia,ja) = iN(2,ia) * jN(1,ja)
-     N3(2,1,2,ia,ja) = iN(1,ia) * jN(2,ja)
-     N3(1,2,2,ia,ja) = iN(1,ia) * jN(2,ja)
-     N3(2,2,2,ia,ja) = iN(0,ia) * jN(3,ja)
-  end do; end do
-  !
-end subroutine TensorBasisFuns
-
-include 'petigarat.f90.in'
-
-end subroutine IGA_BasisFuns_2D
-
-
-pure subroutine IGA_ShapeFuns_2D(&
-     order,                      &
-     nqp,nen,X,                  &
-     M0,M1,M2,M3,                &
-     N0,N1,N2,N3,                &
-     dX,G0,G1,H0,H1)             &
-  bind(C, name="IGA_ShapeFuns_2D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nqp
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: X(dim,nen)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M0(       nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M1(dim,   nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N0(       nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
-  integer(kind=IGA_INTEGER_KIND)  :: q
-  do q=1,nqp
-     call GeometryMap(&
-          order,&
-          nen,X,&
-          M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
-          N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
-  end do
-contains
-include 'petigageo.f90.in'
-end subroutine IGA_ShapeFuns_2D
-
-
-subroutine IGA_BoundaryArea_2D(&
-     m,axis,side,              &
-     geometry,Cx,              &
-     rational,Cw,              &
-     nqp,W,nen,ndr,N,          &
-     dS)                       &
-  bind(C, name="IGA_BoundaryArea_2D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter         :: nsd = 2
-  integer(kind=IGA_INTEGER_KIND), parameter         :: dim = 1
-  integer(kind=IGA_INTEGER_KIND), intent(in)        :: m(2)
-  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: axis, side
-  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: geometry, rational
-  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cx(nsd,m(1),m(2))
-  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cw(    m(1),m(2))
-  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: nqp, nen, ndr
-  real   (kind=IGA_REAL_KIND   ), intent(in)        :: W(nqp), N(0:ndr,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)       :: dS
-  integer(kind=IGA_INTEGER_KIND)  :: k, q
-  real   (kind=IGA_REAL_KIND   )  :: N0(nen), N1(dim,nen), detJ
-  real   (kind=IGA_REAL_KIND   ), pointer :: Xx(:,:), Xw(:)
-  select case (axis)
-  case (0)
-     if (side==0) k=1
-     if (side==1) k=m(1)
-     Xx => Cx(:,k,:); Xw => Cw(k,:)
-  case (1)
-     if (side==0) k=1
-     if (side==1) k=m(2)
-     Xx => Cx(:,:,k); Xw => Cw(:,k)
-  end select
-  detJ = 1
-  dS = 0
-  do q=1,nqp
-     N0(  :) = N(0,:,q)
-     N1(1,:) = N(1,:,q)
-     if (rational /= 0) then
-        call Rationalize(nen,Xw,N0,N1)
-     end if
-     if (geometry /= 0) then
-        call Jacobian(nen,N1,Xx,detJ)
-     end if
-     dS = dS + detJ * W(q)
-  end do
-contains
-pure subroutine Rationalize(nen,W,R0,R1)
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter     :: dim = 1
-  integer(kind=IGA_INTEGER_KIND), intent(in)    :: nen
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: W(nen)
-  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R0(    nen)
-  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R1(dim,nen)
-  integer(kind=IGA_INTEGER_KIND)  :: i
-  real   (kind=IGA_REAL_KIND   )  :: W0
-  R0 = W * R0
-  W0 = sum(R0)
-  R0 = R0 / W0
-  do i=1,dim
-  R1(i,:) = W*R1(i,:) - R0 * sum(W*R1(i,:))
-  end do
-  R1 = R1 / W0
-end subroutine Rationalize
-pure subroutine Jacobian(nen,N,X,J)
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: nsd = 2
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 1
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
-  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
-  real   (kind=IGA_REAL_KIND   ), intent(in)       :: X(nsd,nen)
-  real   (kind=IGA_REAL_KIND   ), intent(out)      :: J
-  real   (kind=IGA_REAL_KIND   )  :: F(dim,nsd), M(dim,dim)
-  F = matmul(N,transpose(X))
-  M = matmul(F,transpose(F))
-  J = sqrt(abs(Determinant(dim,M)))
-end subroutine Jacobian
-include 'petigainv.f90.in'
-end subroutine IGA_BoundaryArea_2D
+pure subroutine IGA_Quadrature_3D(&
+     inq,iX,iW,iJ,                &
+     jnq,jX,jW,jJ,                &
+     knq,kX,kW,kJ,                &
+     X,W,J)                       &
+  bind(C, name="IGA_Quadrature_3D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 3
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jnq
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: knq
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iX(inq), iW(inq), iJ
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jX(jnq), jW(jnq), jJ
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: kX(knq), kW(knq), kJ
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq,jnq,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(    inq,jnq,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(    inq,jnq,knq)
+  integer(kind=IGA_INTEGER_KIND) :: iq
+  integer(kind=IGA_INTEGER_KIND) :: jq
+  integer(kind=IGA_INTEGER_KIND) :: kq
+  do kq=1,knq; do jq=1,jnq; do iq=1,inq
+     X(1,iq,jq,kq) = iX(iq)
+     X(2,iq,jq,kq) = jX(jq)
+     X(3,iq,jq,kq) = kX(kq)
+     W(  iq,jq,kq) = iW(iq) * jW(jq) * kW(kq)
+  end do; end do; end do
+  J = iJ * jJ * kJ
+end subroutine IGA_Quadrature_3D
+
+
+pure subroutine IGA_BasisFuns_3D(&
+     order,                      &
+     rational,W,                 &
+     inq,ina,ind,iN,             &
+     jnq,jna,jnd,jN,             &
+     knq,kna,knd,kN,             &
+     N0,N1,N2,N3)                &
+  bind(C, name="IGA_BasisFuns_3D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 3
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: rational
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq, ina, ind
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jnq, jna, jnd
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: knq, kna, knd
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: W(ina*jna*kna)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: kN(0:knd,kna,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(       ina*jna*kna,inq,jnq,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(   dim,ina*jna*kna,inq,jnq,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(dim**2,ina*jna*kna,inq,jnq,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim**3,ina*jna*kna,inq,jnq,knq)
+  integer(kind=IGA_INTEGER_KIND)  :: ia, iq
+  integer(kind=IGA_INTEGER_KIND)  :: ja, jq
+  integer(kind=IGA_INTEGER_KIND)  :: ka, kq
+  integer(kind=IGA_INTEGER_KIND)  :: nen
+  nen = ina*jna*kna
+  do kq=1,knq
+     do jq=1,jnq
+        do iq=1,inq
+           call TensorBasisFuns(&
+                order,&
+                ina,ind,iN(:,:,iq),&
+                jna,jnd,jN(:,:,jq),&
+                kna,knd,kN(:,:,kq),&
+                N0(  :,iq,jq,kq),&
+                N1(:,:,iq,jq,kq),&
+                N2(:,:,iq,jq,kq),&
+                N3(:,:,iq,jq,kq))
+           if (rational /= 0) then
+              call Rationalize(&
+                   order,&
+                   nen,W,&
+                   N0(  :,iq,jq,kq),&
+                   N1(:,:,iq,jq,kq),&
+                   N2(:,:,iq,jq,kq),&
+                   N3(:,:,iq,jq,kq))
+           end if
+        end do
+     end do
+  end do
+
+contains
+
+pure subroutine TensorBasisFuns(&
+     ord,&
+     ina,ind,iN,&
+     jna,jnd,jN,&
+     kna,knd,kN,&
+     N0,N1,N2,N3)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 3
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ord
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ina, ind
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jna, jnd
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: kna, knd
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: kN(0:knd,kna)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(            ina,jna,kna)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(        dim,ina,jna,kna)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(    dim,dim,ina,jna,kna)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim,dim,dim,ina,jna,kna)
+  integer(kind=IGA_INTEGER_KIND)  :: ia, ja, ka
+  !
+  do ka=1,kna; do ja=1,jna; do ia=1,ina
+     N0(ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(0,ka)
+  end do; end do; end do
+  !
+  do ka=1,kna; do ja=1,jna; do ia=1,ina
+     N1(1,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(0,ka)
+     N1(2,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(0,ka)
+     N1(3,ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(1,ka)
+  end do; end do; end do
+  !
+  if (ord < 2) return ! XXX Optimize!
+  do ka=1,kna; do ja=1,jna; do ia=1,ina
+     N2(1,1,ia,ja,ka) = iN(2,ia) * jN(0,ja) * kN(0,ka)
+     N2(2,1,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(0,ka)
+     N2(3,1,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(1,ka)
+     N2(1,2,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(0,ka)
+     N2(2,2,ia,ja,ka) = iN(0,ia) * jN(2,ja) * kN(0,ka)
+     N2(3,2,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(1,ka)
+     N2(1,3,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(1,ka)
+     N2(2,3,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(1,ka)
+     N2(3,3,ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(2,ka)
+  end do; end do; end do
+  !
+  if (ord < 3) return ! XXX Optimize!
+  do ka=1,kna; do ja=1,jna; do ia=1,ina
+     N3(1,1,1,ia,ja,ka) = iN(3,ia) * jN(0,ja) * kN(0,ka)
+     N3(2,1,1,ia,ja,ka) = iN(2,ia) * jN(1,ja) * kN(0,ka)
+     N3(3,1,1,ia,ja,ka) = iN(2,ia) * jN(0,ja) * kN(1,ka)
+     N3(1,2,1,ia,ja,ka) = iN(2,ia) * jN(1,ja) * kN(0,ka)
+     N3(2,2,1,ia,ja,ka) = iN(1,ia) * jN(2,ja) * kN(0,ka)
+     N3(3,2,1,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
+     N3(1,3,1,ia,ja,ka) = iN(2,ia) * jN(0,ja) * kN(1,ka)
+     N3(2,3,1,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
+     N3(3,3,1,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(2,ka)
+     N3(1,1,2,ia,ja,ka) = iN(2,ia) * jN(1,ja) * kN(0,ka)
+     N3(2,1,2,ia,ja,ka) = iN(1,ia) * jN(2,ja) * kN(0,ka)
+     N3(3,1,2,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
+     N3(1,2,2,ia,ja,ka) = iN(1,ia) * jN(2,ja) * kN(0,ka)
+     N3(2,2,2,ia,ja,ka) = iN(0,ia) * jN(3,ja) * kN(0,ka)
+     N3(3,2,2,ia,ja,ka) = iN(0,ia) * jN(2,ja) * kN(1,ka)
+     N3(1,3,2,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
+     N3(2,3,2,ia,ja,ka) = iN(0,ia) * jN(2,ja) * kN(1,ka)
+     N3(3,3,2,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(2,ka)
+     N3(1,1,3,ia,ja,ka) = iN(2,ia) * jN(0,ja) * kN(1,ka)
+     N3(2,1,3,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
+     N3(3,1,3,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(2,ka)
+     N3(1,2,3,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
+     N3(2,2,3,ia,ja,ka) = iN(0,ia) * jN(2,ja) * kN(1,ka)
+     N3(3,2,3,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(2,ka)
+     N3(1,3,3,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(2,ka)
+     N3(2,3,3,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(2,ka)
+     N3(3,3,3,ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(3,ka)
+  end do; end do; end do
+  !
+end subroutine TensorBasisFuns
+
+include 'petigarat.f90.in'
+
+end subroutine IGA_BasisFuns_3D
+
+
+pure subroutine IGA_ShapeFuns_3D(&
+     order,                      &
+     nqp,nen,X,                  &
+     M0,M1,M2,M3,                &
+     N0,N1,N2,N3,                &
+     dX,G0,G1,H0,H1)             &
+  bind(C, name="IGA_ShapeFuns_3D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 3
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nqp
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: X(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M0(       nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M1(dim,   nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M2(dim**2,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M3(dim**3,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N0(       nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
+  integer(kind=IGA_INTEGER_KIND)  :: q
+  do q=1,nqp
+     call GeometryMap(&
+          order,&
+          nen,X,&
+          M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
+          N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
+          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
+  end do
+contains
+include 'petigageo.f90.in'
+end subroutine IGA_ShapeFuns_3D
+
+
+subroutine IGA_BoundaryArea_3D(&
+     m,axis,side,              &
+     geometry,Cx,              &
+     rational,Cw,              &
+     inqp,iW,inen,indr,iN,     &
+     jnqp,jW,jnen,jndr,jN,     &
+     dS)                       &
+  bind(C, name="IGA_BoundaryArea_3D")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter         :: nsd = 3
+  integer(kind=IGA_INTEGER_KIND), parameter         :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in)        :: m(3)
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: axis, side
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: geometry, rational
+  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cx(nsd,m(1),m(2),m(3))
+  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cw(    m(1),m(2),m(3))
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: inqp, inen, indr
+  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: jnqp, jnen, jndr
+  real   (kind=IGA_REAL_KIND   ), intent(in)        :: iW(inqp), iN(0:indr,inen,inqp)
+  real   (kind=IGA_REAL_KIND   ), intent(in)        :: jW(jnqp), jN(0:jndr,jnen,jnqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)       :: dS
+  integer(kind=IGA_INTEGER_KIND)  :: k, nen, iq, jq, ia, ja
+  real   (kind=IGA_REAL_KIND   )  :: N0(inen,jnen), N1(dim,inen,jnen), detJ
+  real   (kind=IGA_REAL_KIND   ), pointer :: Xx(:,:,:), Xw(:,:)
+  nen = inen*jnen
+  select case (axis)
+  case (0)
+     if (side==0) k=1
+     if (side==1) k=m(1)
+     Xx => Cx(:,k,:,:); Xw => Cw(k,:,:)
+  case (1)
+     if (side==0) k=1
+     if (side==1) k=m(2)
+     Xx => Cx(:,:,k,:); Xw => Cw(:,k,:)
+  case (2)
+     if (side==0) k=1
+     if (side==1) k=m(3)
+     Xx => Cx(:,:,:,k); Xw => Cw(:,:,k)
+  end select
+  detJ = 1
+  dS = 0
+  do jq=1,jnqp
+     do iq=1,inqp
+        do ja=1,jnen; do ia=1,inen
+           N0(ia,ja) = iN(0,ia,iq) * jN(0,ja,jq)
+        end do; end do
+        do ja=1,jnen; do ia=1,inen
+           N1(1,ia,ja) = iN(1,ia,iq) * jN(0,ja,jq)
+           N1(2,ia,ja) = iN(0,ia,iq) * jN(1,ja,jq)
+        end do; end do
+        if (rational /= 0) then
+           call Rationalize(nen,Xw,N0,N1)
+        end if
+        if (geometry /= 0) then
+           call Jacobian(nen,N1,Xx,detJ)
+        end if
+        dS = dS + detJ * (iW(iq)*jW(jq))
+     end do
+  end do
+contains
+pure subroutine Rationalize(nen,W,R0,R1)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter     :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in)    :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)    :: W(nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R0(    nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R1(dim,nen)
+  integer(kind=IGA_INTEGER_KIND)  :: i
+  real   (kind=IGA_REAL_KIND   )  :: W0
+  R0 = W * R0
+  W0 = sum(R0)
+  R0 = R0 / W0
+  do i=1,dim
+  R1(i,:) = W*R1(i,:) - R0 * sum(W*R1(i,:))
+  end do
+  R1 = R1 / W0
+end subroutine Rationalize
+pure subroutine Jacobian(nen,N,X,J)
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), parameter        :: nsd = 3
+  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: X(nsd,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: J
+  real   (kind=IGA_REAL_KIND   )  :: F(dim,nsd), M(dim,dim)
+  F = matmul(N,transpose(X))
+  M = matmul(F,transpose(F))
+  J = sqrt(abs(Determinant(dim,M)))
+end subroutine Jacobian
+include 'petigainv.f90.in'
+end subroutine IGA_BoundaryArea_3D

src/petiga3d.f90

-pure subroutine IGA_Quadrature_3D(&
-     inq,iX,iW,iJ,                &
-     jnq,jX,jW,jJ,                &
-     knq,kX,kW,kJ,                &
-     X,W,J)                       &
-  bind(C, name="IGA_Quadrature_3D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 3
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jnq
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: knq
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iX(inq), iW(inq), iJ
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jX(jnq), jW(jnq), jJ
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: kX(knq), kW(knq), kJ
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq,jnq,knq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(    inq,jnq,knq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(    inq,jnq,knq)
-  integer(kind=IGA_INTEGER_KIND) :: iq
-  integer(kind=IGA_INTEGER_KIND) :: jq
-  integer(kind=IGA_INTEGER_KIND) :: kq
-  do kq=1,knq; do jq=1,jnq; do iq=1,inq
-     X(1,iq,jq,kq) = iX(iq)
-     X(2,iq,jq,kq) = jX(jq)
-     X(3,iq,jq,kq) = kX(kq)
-     W(  iq,jq,kq) = iW(iq) * jW(jq) * kW(kq)
-  end do; end do; end do
-  J = iJ * jJ * kJ
-end subroutine IGA_Quadrature_3D
-
-
-pure subroutine IGA_BasisFuns_3D(&
-     order,                      &
-     rational,W,                 &
-     inq,ina,ind,iN,             &
-     jnq,jna,jnd,jN,             &
-     knq,kna,knd,kN,             &
-     N0,N1,N2,N3)                &
-  bind(C, name="IGA_BasisFuns_3D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 3
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: rational
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: inq, ina, ind
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jnq, jna, jnd
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: knq, kna, knd
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: W(ina*jna*kna)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina,inq)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna,jnq)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: kN(0:knd,kna,knq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(       ina*jna*kna,inq,jnq,knq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(   dim,ina*jna*kna,inq,jnq,knq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(dim**2,ina*jna*kna,inq,jnq,knq)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim**3,ina*jna*kna,inq,jnq,knq)
-  integer(kind=IGA_INTEGER_KIND)  :: ia, iq
-  integer(kind=IGA_INTEGER_KIND)  :: ja, jq
-  integer(kind=IGA_INTEGER_KIND)  :: ka, kq
-  integer(kind=IGA_INTEGER_KIND)  :: nen
-  nen = ina*jna*kna
-  do kq=1,knq
-     do jq=1,jnq
-        do iq=1,inq
-           call TensorBasisFuns(&
-                order,&
-                ina,ind,iN(:,:,iq),&
-                jna,jnd,jN(:,:,jq),&
-                kna,knd,kN(:,:,kq),&
-                N0(  :,iq,jq,kq),&
-                N1(:,:,iq,jq,kq),&
-                N2(:,:,iq,jq,kq),&
-                N3(:,:,iq,jq,kq))
-           if (rational /= 0) then
-              call Rationalize(&
-                   order,&
-                   nen,W,&
-                   N0(  :,iq,jq,kq),&
-                   N1(:,:,iq,jq,kq),&
-                   N2(:,:,iq,jq,kq),&
-                   N3(:,:,iq,jq,kq))
-           end if
-        end do
-     end do
-  end do
-
-contains
-
-pure subroutine TensorBasisFuns(&
-     ord,&
-     ina,ind,iN,&
-     jna,jnd,jN,&
-     kna,knd,kN,&
-     N0,N1,N2,N3)
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 3
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ord
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: ina, ind
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: jna, jnd
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: kna, knd
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: iN(0:ind,ina)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna)
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: kN(0:knd,kna)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N0(            ina,jna,kna)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N1(        dim,ina,jna,kna)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N2(    dim,dim,ina,jna,kna)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: N3(dim,dim,dim,ina,jna,kna)
-  integer(kind=IGA_INTEGER_KIND)  :: ia, ja, ka
-  !
-  do ka=1,kna; do ja=1,jna; do ia=1,ina
-     N0(ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(0,ka)
-  end do; end do; end do
-  !
-  do ka=1,kna; do ja=1,jna; do ia=1,ina
-     N1(1,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(0,ka)
-     N1(2,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(0,ka)
-     N1(3,ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(1,ka)
-  end do; end do; end do
-  !
-  if (ord < 2) return ! XXX Optimize!
-  do ka=1,kna; do ja=1,jna; do ia=1,ina
-     N2(1,1,ia,ja,ka) = iN(2,ia) * jN(0,ja) * kN(0,ka)
-     N2(2,1,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(0,ka)
-     N2(3,1,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(1,ka)
-     N2(1,2,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(0,ka)
-     N2(2,2,ia,ja,ka) = iN(0,ia) * jN(2,ja) * kN(0,ka)
-     N2(3,2,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(1,ka)
-     N2(1,3,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(1,ka)
-     N2(2,3,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(1,ka)
-     N2(3,3,ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(2,ka)
-  end do; end do; end do
-  !
-  if (ord < 3) return ! XXX Optimize!
-  do ka=1,kna; do ja=1,jna; do ia=1,ina
-     N3(1,1,1,ia,ja,ka) = iN(3,ia) * jN(0,ja) * kN(0,ka)
-     N3(2,1,1,ia,ja,ka) = iN(2,ia) * jN(1,ja) * kN(0,ka)
-     N3(3,1,1,ia,ja,ka) = iN(2,ia) * jN(0,ja) * kN(1,ka)
-     N3(1,2,1,ia,ja,ka) = iN(2,ia) * jN(1,ja) * kN(0,ka)
-     N3(2,2,1,ia,ja,ka) = iN(1,ia) * jN(2,ja) * kN(0,ka)
-     N3(3,2,1,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
-     N3(1,3,1,ia,ja,ka) = iN(2,ia) * jN(0,ja) * kN(1,ka)
-     N3(2,3,1,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
-     N3(3,3,1,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(2,ka)
-     N3(1,1,2,ia,ja,ka) = iN(2,ia) * jN(1,ja) * kN(0,ka)
-     N3(2,1,2,ia,ja,ka) = iN(1,ia) * jN(2,ja) * kN(0,ka)
-     N3(3,1,2,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
-     N3(1,2,2,ia,ja,ka) = iN(1,ia) * jN(2,ja) * kN(0,ka)
-     N3(2,2,2,ia,ja,ka) = iN(0,ia) * jN(3,ja) * kN(0,ka)
-     N3(3,2,2,ia,ja,ka) = iN(0,ia) * jN(2,ja) * kN(1,ka)
-     N3(1,3,2,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
-     N3(2,3,2,ia,ja,ka) = iN(0,ia) * jN(2,ja) * kN(1,ka)
-     N3(3,3,2,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(2,ka)
-     N3(1,1,3,ia,ja,ka) = iN(2,ia) * jN(0,ja) * kN(1,ka)
-     N3(2,1,3,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
-     N3(3,1,3,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(2,ka)
-     N3(1,2,3,ia,ja,ka) = iN(1,ia) * jN(1,ja) * kN(1,ka)
-     N3(2,2,3,ia,ja,ka) = iN(0,ia) * jN(2,ja) * kN(1,ka)
-     N3(3,2,3,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(2,ka)
-     N3(1,3,3,ia,ja,ka) = iN(1,ia) * jN(0,ja) * kN(2,ka)
-     N3(2,3,3,ia,ja,ka) = iN(0,ia) * jN(1,ja) * kN(2,ka)
-     N3(3,3,3,ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(3,ka)
-  end do; end do; end do
-  !
-end subroutine TensorBasisFuns
-
-include 'petigarat.f90.in'
-
-end subroutine IGA_BasisFuns_3D
-
-
-pure subroutine IGA_ShapeFuns_3D(&
-     order,                      &
-     nqp,nen,X,                  &
-     M0,M1,M2,M3,                &
-     N0,N1,N2,N3,                &
-     dX,G0,G1,H0,H1)             &
-  bind(C, name="IGA_ShapeFuns_3D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 3
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: order
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nqp
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: X(dim,nen)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M0(       nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M1(dim,   nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: M3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N0(       nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N1(dim,   nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: dX(nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G0(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G1(dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H0(dim,dim,dim,nqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)   :: H1(dim,dim,dim,nqp)
-  integer(kind=IGA_INTEGER_KIND)  :: q
-  do q=1,nqp
-     call GeometryMap(&
-          order,&
-          nen,X,&
-          M0(:,q),M1(:,:,q),M2(:,:,q),M3(:,:,q),&
-          N0(:,q),N1(:,:,q),N2(:,:,q),N3(:,:,q),&
-          dX(q),G0(:,:,q),G1(:,:,q),H0(:,:,:,q),H1(:,:,:,q))
-  end do
-contains
-include 'petigageo.f90.in'
-end subroutine IGA_ShapeFuns_3D
-
-
-subroutine IGA_BoundaryArea_3D(&
-     m,axis,side,              &
-     geometry,Cx,              &
-     rational,Cw,              &
-     inqp,iW,inen,indr,iN,     &
-     jnqp,jW,jnen,jndr,jN,     &
-     dS)                       &
-  bind(C, name="IGA_BoundaryArea_3D")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter         :: nsd = 3
-  integer(kind=IGA_INTEGER_KIND), parameter         :: dim = 2
-  integer(kind=IGA_INTEGER_KIND), intent(in)        :: m(3)
-  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: axis, side
-  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: geometry, rational
-  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cx(nsd,m(1),m(2),m(3))
-  real   (kind=IGA_REAL_KIND   ), intent(in),target :: Cw(    m(1),m(2),m(3))
-  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: inqp, inen, indr
-  integer(kind=IGA_INTEGER_KIND), intent(in),value  :: jnqp, jnen, jndr
-  real   (kind=IGA_REAL_KIND   ), intent(in)        :: iW(inqp), iN(0:indr,inen,inqp)
-  real   (kind=IGA_REAL_KIND   ), intent(in)        :: jW(jnqp), jN(0:jndr,jnen,jnqp)
-  real   (kind=IGA_REAL_KIND   ), intent(out)       :: dS
-  integer(kind=IGA_INTEGER_KIND)  :: k, nen, iq, jq, ia, ja
-  real   (kind=IGA_REAL_KIND   )  :: N0(inen,jnen), N1(dim,inen,jnen), detJ
-  real   (kind=IGA_REAL_KIND   ), pointer :: Xx(:,:,:), Xw(:,:)
-  nen = inen*jnen
-  select case (axis)
-  case (0)
-     if (side==0) k=1
-     if (side==1) k=m(1)
-     Xx => Cx(:,k,:,:); Xw => Cw(k,:,:)
-  case (1)
-     if (side==0) k=1
-     if (side==1) k=m(2)
-     Xx => Cx(:,:,k,:); Xw => Cw(:,k,:)
-  case (2)
-     if (side==0) k=1
-     if (side==1) k=m(3)
-     Xx => Cx(:,:,:,k); Xw => Cw(:,:,k)
-  end select
-  detJ = 1
-  dS = 0
-  do jq=1,jnqp
-     do iq=1,inqp
-        do ja=1,jnen; do ia=1,inen
-           N0(ia,ja) = iN(0,ia,iq) * jN(0,ja,jq)
-        end do; end do
-        do ja=1,jnen; do ia=1,inen
-           N1(1,ia,ja) = iN(1,ia,iq) * jN(0,ja,jq)
-           N1(2,ia,ja) = iN(0,ia,iq) * jN(1,ja,jq)
-        end do; end do
-        if (rational /= 0) then
-           call Rationalize(nen,Xw,N0,N1)
-        end if
-        if (geometry /= 0) then
-           call Jacobian(nen,N1,Xx,detJ)
-        end if
-        dS = dS + detJ * (iW(iq)*jW(jq))
-     end do
-  end do
-contains
-pure subroutine Rationalize(nen,W,R0,R1)
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter     :: dim = 2
-  integer(kind=IGA_INTEGER_KIND), intent(in)    :: nen
-  real   (kind=IGA_REAL_KIND   ), intent(in)    :: W(nen)
-  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R0(    nen)
-  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R1(dim,nen)
-  integer(kind=IGA_INTEGER_KIND)  :: i
-  real   (kind=IGA_REAL_KIND   )  :: W0
-  R0 = W * R0
-  W0 = sum(R0)
-  R0 = R0 / W0
-  do i=1,dim
-  R1(i,:) = W*R1(i,:) - R0 * sum(W*R1(i,:))
-  end do
-  R1 = R1 / W0
-end subroutine Rationalize
-pure subroutine Jacobian(nen,N,X,J)
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), parameter        :: nsd = 3
-  integer(kind=IGA_INTEGER_KIND), parameter        :: dim = 2
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen
-  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
-  real   (kind=IGA_REAL_KIND   ), intent(in)       :: X(nsd,nen)
-  real   (kind=IGA_REAL_KIND   ), intent(out)      :: J
-  real   (kind=IGA_REAL_KIND   )  :: F(dim,nsd), M(dim,dim)
-  F = matmul(N,transpose(X))
-  M = matmul(F,transpose(F))
-  J = sqrt(abs(Determinant(dim,M)))
-end subroutine Jacobian
-include 'petigainv.f90.in'
-end subroutine IGA_BoundaryArea_3D

src/petigabsp.F90

+! -*- f90 -*-
+
+subroutine IGA_Basis_BSpline(k,uu,p,d,U,B) &
+  bind(C, name="IGA_Basis_BSpline")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: k, p, d
+  real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:k+p)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
+  real   (kind=IGA_REAL_KIND   )  :: ders(0:p,0:d)
+  call BasisFunsDers(k,uu,p,d,U,ders)
+  B = transpose(ders)
+contains
+pure subroutine BasisFunsDers(i,uu,p,n,U,ders)
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), intent(in)  :: i, p, n
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: uu, U(0:i+p)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: ders(0:p,0:n)
+  integer(kind=IGA_INTEGER_KIND)  :: j, k, r, s1, s2, rk, pk, j1, j2
+  real   (kind=IGA_REAL_KIND   )  :: saved, temp, d
+  real   (kind=IGA_REAL_KIND   )  :: left(p), right(p)
+  real   (kind=IGA_REAL_KIND   )  :: ndu(0:p,0:p), a(0:1,0:p)
+  ndu(0,0) = 1
+  do j = 1, p
+     left(j)  = uu - U(i+1-j)
+     right(j) = U(i+j) - uu
+     saved = 0
+     do r = 0, j-1
+        ndu(j,r) = right(r+1) + left(j-r)
+        temp = ndu(r,j-1) / ndu(j,r)
+        ndu(r,j) = saved + right(r+1) * temp
+        saved = left(j-r) * temp
+     end do
+     ndu(j,j) = saved
+  end do
+  ders(:,0) = ndu(:,p)
+  do r = 0, p
+     s1 = 0; s2 = 1;
+     a(0,0) = 1
+     do k = 1, n
+        d = 0
+        rk = r-k; pk = p-k;
+        if (r >= k) then
+           a(s2,0) = a(s1,0) / ndu(pk+1,rk)
+           d =  a(s2,0) * ndu(rk,pk)
+        end if
+        if (rk > -1) then
+           j1 = 1
+        else
+           j1 = -rk
+        end if
+        if (r-1 <= pk) then
+           j2 = k-1
+        else
+           j2 = p-r
+        end if
+        do j = j1, j2
+           a(s2,j) = (a(s1,j) - a(s1,j-1)) / ndu(pk+1,rk+j)
+           d =  d + a(s2,j) * ndu(rk+j,pk)
+        end do
+        if (r <= pk) then
+           a(s2,k) = - a(s1,k-1) / ndu(pk+1,r)
+           d =  d + a(s2,k) * ndu(r,pk)
+        end if
+        ders(r,k) = d
+        j = s1; s1 = s2; s2 = j;
+     end do
+  end do
+  r = p
+  do k = 1, n
+     ders(:,k) = ders(:,k) * r
+     r = r * (p-k)
+  end do
+end subroutine BasisFunsDers
+end subroutine IGA_Basis_BSpline
+
+
+subroutine IGA_Basis_Lagrange(kk,uu,p,d,U,B) &
+  bind(C, name="IGA_Basis_Lagrange")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: kk, p, d
+  real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:kk+p)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
+  integer(kind=IGA_INTEGER_KIND)  :: m, i, j, k, l
+  real   (kind=IGA_REAL_KIND   )  :: Lp, Ls1, Ls2, Ls3
+  real   (kind=IGA_REAL_KIND   )  :: X(0:p)
+
+  do m = 0, p
+     X(m) = U(kk) + m * (U(kk+1) - U(kk)) / p
+  end do
+
+  do m = 0, p
+     Lp = 1
+     do i = 0, p
+        if (i == m) cycle
+        Lp = Lp * (uu-X(i))/(X(m)-X(i))
+     end do
+     B(0,m) = Lp
+  end do
+
+  if (d < 1) return
+  do m = 0, p
+     Ls1 = 0
+     do j = 0, p
+        if (j == m) cycle
+        Lp = 1
+        do i = 0, p
+           if (i == m .or. i == j) cycle
+           Lp = Lp * (uu-X(i))/(X(m)-X(i))
+        end do
+        Ls1 = Ls1 + Lp/(X(m)-X(j))
+     end do
+     B(1,m) = Ls1
+  end do
+
+  if (d < 2) return
+  do m = 0, p
+     Ls2 = 0
+     do k = 0, p
+        if (k == m) cycle
+        Ls1 = 0
+        do j = 0, p
+           if (j == m .or. j == k) cycle
+           Lp = 1
+           do i = 0, p
+              if (i == m .or. i == k .or. i == j) cycle
+              Lp = Lp * (uu-X(i))/(X(m)-X(i))
+           end do
+           Ls1 = Ls1 + Lp/(X(m)-X(j))
+        end do
+        Ls2 = Ls2 + Ls1/(X(m)-X(k))
+     end do
+     B(2,m) = Ls2
+  end do
+
+  if (d < 3) return
+  do m = 0, p
+     Ls3 = 0
+     do l = 0, p
+        if (l == m) cycle
+        Ls2 = 0
+        do k = 0, p
+           if (k == m .or. k == l) cycle
+           Ls1 = 0
+           do j = 0, p
+              if (j == m .or. j == l .or. j == k) cycle
+              Lp = 1
+              do i = 0, p
+                 if (i == m .or. i == l .or. i == k .or. i == j) cycle
+                 Lp = Lp * (uu-X(i))/(X(m)-X(i))
+              end do
+              Ls1 = Ls1 + Lp/(X(m)-X(j))
+           end do
+           Ls2 = Ls2 + Ls1/(X(m)-X(k))
+        end do
+        Ls3 = Ls3 + Ls2/(X(m)-X(l))
+     end do
+     B(3,m) = Ls3
+  end do
+
+end subroutine IGA_Basis_Lagrange
+
+
+subroutine IGA_Basis_Hierarchical(kk,uu,p,d,U,B) &
+  bind(C, name="IGA_Basis_Hierarchical")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: kk, p, d
+  real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:kk+p)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
+  integer(kind=IGA_INTEGER_KIND)  :: i, k
+  real   (kind=IGA_REAL_KIND   )  :: J, x, Lp(0:p,0:d)
+  real   (kind=IGA_REAL_KIND   ), parameter :: two = 2
+
+  J = (U(kk+1)-U(kk))/2
+  x = (uu-U(kk))/J - 1
+
+  B(0,0) = (1-x)/2
+  B(0,p) = (x+1)/2
+  if (d > 0) then
+     B(1,0) = -1/two
+     B(1,p) = +1/two
+  endif
+
+  if (p > 1) then
+     Lp(:,:) = 0
+     Lp(0,0) = 1
+     Lp(1,0) = x
+     if (d > 0) then
+        Lp(0,1) = 0
+        Lp(1,1) = 1
+     end if
+     do i = 1, p-1
+        Lp(i+1,0) = ((2*i+1)*x*Lp(i,0) - i*Lp(i-1,0))/(i+1)
+        B(0,i) = (-Lp(i+1,0) + Lp(i-1,0))/sqrt(two*(2*(i+1)-1))
+        do k = 1, d
+           Lp(i+1,k) = (2*i+1)*Lp(i,k-1) + Lp(i-1,k)
+           B(k,i) = (-Lp(i+1,k) + Lp(i-1,k))/sqrt(two*(2*(i+1)-1))
+        end do
+     end do
+  end if
+
+  do k = 1, d
+     B(k,:) = B(k,:)/(J**k)
+  end do
+
+end subroutine IGA_Basis_Hierarchical

src/petigabsp.f90

-! -*- f90 -*-
-
-subroutine IGA_Basis_BSpline(k,uu,p,d,U,B) &
-  bind(C, name="IGA_Basis_BSpline")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: k, p, d
-  real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
-  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:k+p)
-  real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
-  real   (kind=IGA_REAL_KIND   )  :: ders(0:p,0:d)
-  call BasisFunsDers(k,uu,p,d,U,ders)
-  B = transpose(ders)
-contains
-pure subroutine BasisFunsDers(i,uu,p,n,U,ders)
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), intent(in)  :: i, p, n
-  real   (kind=IGA_REAL_KIND   ), intent(in)  :: uu, U(0:i+p)
-  real   (kind=IGA_REAL_KIND   ), intent(out) :: ders(0:p,0:n)
-  integer(kind=IGA_INTEGER_KIND)  :: j, k, r, s1, s2, rk, pk, j1, j2
-  real   (kind=IGA_REAL_KIND   )  :: saved, temp, d
-  real   (kind=IGA_REAL_KIND   )  :: left(p), right(p)
-  real   (kind=IGA_REAL_KIND   )  :: ndu(0:p,0:p), a(0:1,0:p)
-  ndu(0,0) = 1
-  do j = 1, p
-     left(j)  = uu - U(i+1-j)
-     right(j) = U(i+j) - uu
-     saved = 0
-     do r = 0, j-1
-        ndu(j,r) = right(r+1) + left(j-r)
-        temp = ndu(r,j-1) / ndu(j,r)
-        ndu(r,j) = saved + right(r+1) * temp
-        saved = left(j-r) * temp
-     end do
-     ndu(j,j) = saved
-  end do
-  ders(:,0) = ndu(:,p)
-  do r = 0, p
-     s1 = 0; s2 = 1;
-     a(0,0) = 1
-     do k = 1, n
-        d = 0
-        rk = r-k; pk = p-k;
-        if (r >= k) then
-           a(s2,0) = a(s1,0) / ndu(pk+1,rk)
-           d =  a(s2,0) * ndu(rk,pk)
-        end if
-        if (rk > -1) then
-           j1 = 1
-        else
-           j1 = -rk
-        end if
-        if (r-1 <= pk) then
-           j2 = k-1
-        else
-           j2 = p-r
-        end if
-        do j = j1, j2
-           a(s2,j) = (a(s1,j) - a(s1,j-1)) / ndu(pk+1,rk+j)
-           d =  d + a(s2,j) * ndu(rk+j,pk)
-        end do
-        if (r <= pk) then
-           a(s2,k) = - a(s1,k-1) / ndu(pk+1,r)
-           d =  d + a(s2,k) * ndu(r,pk)
-        end if
-        ders(r,k) = d
-        j = s1; s1 = s2; s2 = j;
-     end do
-  end do
-  r = p
-  do k = 1, n
-     ders(:,k) = ders(:,k) * r
-     r = r * (p-k)
-  end do
-end subroutine BasisFunsDers
-end subroutine IGA_Basis_BSpline
-
-
-subroutine IGA_Basis_Lagrange(kk,uu,p,d,U,B) &
-  bind(C, name="IGA_Basis_Lagrange")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: kk, p, d
-  real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
-  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:kk+p)
-  real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
-  integer(kind=IGA_INTEGER_KIND)  :: m, i, j, k, l
-  real   (kind=IGA_REAL_KIND   )  :: Lp, Ls1, Ls2, Ls3
-  real   (kind=IGA_REAL_KIND   )  :: X(0:p)
-
-  do m = 0, p
-     X(m) = U(kk) + m * (U(kk+1) - U(kk)) / p
-  end do
-
-  do m = 0, p
-     Lp = 1
-     do i = 0, p
-        if (i == m) cycle
-        Lp = Lp * (uu-X(i))/(X(m)-X(i))
-     end do
-     B(0,m) = Lp
-  end do
-
-  if (d < 1) return
-  do m = 0, p
-     Ls1 = 0
-     do j = 0, p
-        if (j == m) cycle
-        Lp = 1
-        do i = 0, p
-           if (i == m .or. i == j) cycle
-           Lp = Lp * (uu-X(i))/(X(m)-X(i))
-        end do
-        Ls1 = Ls1 + Lp/(X(m)-X(j))
-     end do
-     B(1,m) = Ls1
-  end do
-
-  if (d < 2) return
-  do m = 0, p
-     Ls2 = 0
-     do k = 0, p
-        if (k == m) cycle
-        Ls1 = 0
-        do j = 0, p
-           if (j == m .or. j == k) cycle
-           Lp = 1
-           do i = 0, p
-              if (i == m .or. i == k .or. i == j) cycle
-              Lp = Lp * (uu-X(i))/(X(m)-X(i))
-           end do
-           Ls1 = Ls1 + Lp/(X(m)-X(j))
-        end do
-        Ls2 = Ls2 + Ls1/(X(m)-X(k))
-     end do
-     B(2,m) = Ls2
-  end do
-
-  if (d < 3) return
-  do m = 0, p
-     Ls3 = 0
-     do l = 0, p
-        if (l == m) cycle
-        Ls2 = 0
-        do k = 0, p
-           if (k == m .or. k == l) cycle
-           Ls1 = 0
-           do j = 0, p
-              if (j == m .or. j == l .or. j == k) cycle
-              Lp = 1
-              do i = 0, p
-                 if (i == m .or. i == l .or. i == k .or. i == j) cycle
-                 Lp = Lp * (uu-X(i))/(X(m)-X(i))
-              end do
-              Ls1 = Ls1 + Lp/(X(m)-X(j))
-           end do
-           Ls2 = Ls2 + Ls1/(X(m)-X(k))
-        end do
-        Ls3 = Ls3 + Ls2/(X(m)-X(l))
-     end do
-     B(3,m) = Ls3
-  end do
-
-end subroutine IGA_Basis_Lagrange
-
-
-subroutine IGA_Basis_Hierarchical(kk,uu,p,d,U,B) &
-  bind(C, name="IGA_Basis_Hierarchical")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INTEGER_KIND), intent(in),value :: kk, p, d
-  real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
-  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:kk+p)
-  real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
-  integer(kind=IGA_INTEGER_KIND)  :: i, k
-  real   (kind=IGA_REAL_KIND   )  :: J, x, Lp(0:p,0:d)
-  real   (kind=IGA_REAL_KIND   ), parameter :: two = 2
-
-  J = (U(kk+1)-U(kk))/2
-  x = (uu-U(kk))/J - 1
-
-  B(0,0) = (1-x)/2
-  B(0,p) = (x+1)/2
-  if (d > 0) then
-     B(1,0) = -1/two
-     B(1,p) = +1/two
-  endif
-
-  if (p > 1) then
-     Lp(:,:) = 0
-     Lp(0,0) = 1
-     Lp(1,0) = x
-     if (d > 0) then
-        Lp(0,1) = 0
-        Lp(1,1) = 1
-     end if
-     do i = 1, p-1
-        Lp(i+1,0) = ((2*i+1)*x*Lp(i,0) - i*Lp(i-1,0))/(i+1)
-        B(0,i) = (-Lp(i+1,0) + Lp(i-1,0))/sqrt(two*(2*(i+1)-1))
-        do k = 1, d
-           Lp(i+1,k) = (2*i+1)*Lp(i,k-1) + Lp(i-1,k)
-           B(k,i) = (-Lp(i+1,k) + Lp(i-1,k))/sqrt(two*(2*(i+1)-1))
-        end do
-     end do
-  end if
-
-  do k = 1, d
-     B(k,:) = B(k,:)/(J**k)
-  end do
-
-end subroutine IGA_Basis_Hierarchical

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)
+  integer, parameter :: rk = IGA_REAL_KIND
+  select case (q)
+  case (1) ! p <= 1
+     X(0) = 0.0_rk
+     W(0) = 2.0_rk
+  case (2) ! p <= 3
+     X(0) = -0.577350269189625764509148780501957455_rk ! 1/sqrt(3)
+     X(1) = -X(0)
+     W(0) =  1.0_rk                                    ! 1
+     W(1) =  W(0)
+  case (3) ! p <= 5
+     X(0) = -0.774596669241483377035853079956479922_rk ! sqrt(3/5)
+     X(1) =  0.0_rk                                    ! 0
+     X(2) = -X(0)
+     W(0) =  0.555555555555555555555555555555555556_rk ! 5/9
+     W(1) =  0.888888888888888888888888888888888889_rk ! 8/9
+     W(2) =  W(0)
+  case (4) ! p <= 7
+     X(0) = -0.861136311594052575223946488892809506_rk ! sqrt((3+2*sqrt(6/5))/7)
+     X(1) = -0.339981043584856264802665759103244686_rk ! sqrt((3-2*sqrt(6/5))/7)
+     X(2) = -X(1)
+     X(3) = -X(0)
+     W(0) =  0.347854845137453857373063949221999408_rk ! (18-sqrt(30))/36
+     W(1) =  0.652145154862546142626936050778000592_rk ! (18+sqrt(30))/36
+     W(2) =  W(1)
+     W(3) =  W(0)
+  case (5) ! p <= 9
+     X(0) = -0.906179845938663992797626878299392962_rk ! 1/3*sqrt(5+2*sqrt(10/7))
+     X(1) = -0.538469310105683091036314420700208806_rk ! 1/3*sqrt(5-2*sqrt(10/7))
+     X(2) =  0.0_rk                                    ! 0
+     X(3) = -X(1)
+     X(4) = -X(0)
+     W(0) =  0.236926885056189087514264040719917362_rk ! (322-13*sqrt(70))/900
+     W(1) =  0.478628670499366468041291514835638193_rk ! (322+13*sqrt(70))/900
+     W(2) =  0.568888888888888888888888888888888889_rk ! 128/225
+     W(3) =  W(1)
+     W(4) =  W(0)
+  case (6) ! p <= 11
+     X(0) = -0.9324695142031520278123015544939946_rk ! << NumericalDifferentialEquationAnalysis`
+     X(1) = -0.6612093864662645136613995950199053_rk ! GaussianQuadratureWeights(6, -1, 1, 37)
+     X(2) = -0.2386191860831969086305017216807119_rk
+     X(3) = -X(2)
+     X(4) = -X(1)
+     X(5) = -X(0)
+     W(0) =  0.171324492379170345040296142172732894_rk
+     W(1) =  0.360761573048138607569833513837716112_rk
+     W(2) =  0.467913934572691047389870343989550995_rk
+     W(3) =  W(2)
+     W(4) =  W(1)
+     W(5) =  W(0)
+  case (7) ! p <= 13
+     X(0) = -0.9491079123427585245261896840478513_rk ! << NumericalDifferentialEquationAnalysis`
+     X(1) = -0.7415311855993944398638647732807884_rk ! GaussianQuadratureWeights(7, -1, 1, 37)
+     X(2) = -0.4058451513773971669066064120769615_rk
+     X(3) =  0.0_rk
+     X(4) = -X(2)
+     X(5) = -X(1)
+     X(6) = -X(0)
+     W(0) =  0.129484966168869693270611432679082018_rk
+     W(1) =  0.279705391489276667901467771423779582_rk
+     W(2) =  0.381830050505118944950369775488975134_rk
+     W(3) =  0.417959183673469387755102040816326531_rk
+     W(4) =  W(2)
+     W(5) =  W(1)
+     W(6) =  W(0)
+  case (8) ! p <= 15
+     X(0) = -0.9602898564975362316835608685694730_rk ! << NumericalDifferentialEquationAnalysis`
+     X(1) = -0.7966664774136267395915539364758304_rk ! GaussianQuadratureWeights(8, -1, 1, 37)
+     X(2) = -0.5255324099163289858177390491892463_rk
+     X(3) = -0.1834346424956498049394761423601840_rk
+     X(4) = -X(3)
+     X(5) = -X(2)
+     X(6) = -X(1)
+     X(7) = -X(0)
+     W(0) =  0.101228536290376259152531354309962190_rk
+     W(1) =  0.222381034453374470544355994426240884_rk
+     W(2) =  0.313706645877887287337962201986601313_rk
+     W(3) =  0.362683783378361982965150449277195612_rk
+     W(4) =  W(3)
+     W(5) =  W(2)
+     W(6) =  W(1)
+     W(7) =  W(0)
+  case (9) ! p <= 17
+     X(0) = -0.9681602395076260898355762029036729_rk ! << NumericalDifferentialEquationAnalysis`
+     X(1) = -0.8360311073266357942994297880697349_rk ! GaussianQuadratureWeights(9, -1, 1, 37)
+     X(2) = -0.6133714327005903973087020393414742_rk
+     X(3) = -0.3242534234038089290385380146433366_rk
+     X(4) =  0.0_rk
+     X(5) = -X(3)
+     X(6) = -X(2)
+     X(7) = -X(1)
+     X(8) = -X(0)
+     W(0) =  0.081274388361574411971892158110523651_rk
+     W(1) =  0.180648160694857404058472031242912810_rk
+     W(2) =  0.260610696402935462318742869418632850_rk
+     W(3) =  0.312347077040002840068630406584443666_rk
+     W(4) =  0.330239355001259763164525069286974049_rk
+     W(5) =  W(3)
+     W(6) =  W(2)
+     W(7) =  W(1)
+     W(8) =  W(0)
+  case default
+     X = 0.0_rk
+     W = 0.0_rk
+  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)
+  integer, parameter :: rk = IGA_REAL_KIND
+  select case (q)
+  case (2) ! p <= 1
+     X(0) = -1.0_rk
+     X(1) = -X(0)
+     W(0) =  1.0_rk
+     W(1) =  W(0)
+  case (3) ! p <= 3
+     X(0) = -1.0_rk                                    ! -1
+     X(1) =  0.0_rk                                    ! 0
+     X(2) = -X(0)
+     W(0) =  0.333333333333333333333333333333333333_rk ! 1/3
+     W(1) =  1.333333333333333333333333333333333333_rk ! 4/3
+     W(2) =  W(0)
+  case (4) ! p <= 5
+     X(0) = -1.0_rk                                    ! -1
+     X(1) = -0.447213595499957939281834733746255246_rk ! 1/sqrt(5)
+     X(2) = -X(1)
+     X(3) = -X(0)
+     W(0) =  0.166666666666666666666666666666666667_rk ! 1/6
+     W(1) =  0.833333333333333333333333333333333343_rk ! 5/6
+     W(2) =  W(1)
+     W(3) =  W(0)
+  case (5) ! p <= 7
+     X(0) = -1.0_rk                                    ! -1
+     X(1) = -0.654653670707977143798292456246858356_rk ! sqrt(3/7)
+     X(2) =  0.0_rk                                    ! 0
+     X(3) = -X(1)
+     X(4) = -X(0)
+     W(0) =  0.1_rk                                    !  1/10
+     W(1) =  0.544444444444444444444444444444444444_rk ! 49/90
+     W(2) =  0.711111111111111111111111111111111111_rk ! 32/45
+     W(3) =  W(1)
+     W(4) =  W(0)
+  case default
+     X = 0.0_rk
+     W = 0.0_rk
+  end select
+end subroutine IGA_Rule_GaussLobatto

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)
-  integer, parameter :: rk = IGA_REAL_KIND
-  select case (q)
-  case (1) ! p <= 1
-     X(0) = 0.0_rk
-     W(0) = 2.0_rk
-  case (2) ! p <= 3
-     X(0) = -0.577350269189625764509148780501957455_rk ! 1/sqrt(3)
-     X(1) = -X(0)
-     W(0) =  1.0_rk                                    ! 1
-     W(1) =  W(0)
-  case (3) ! p <= 5
-     X(0) = -0.774596669241483377035853079956479922_rk ! sqrt(3/5)
-     X(1) =  0.0_rk                                    ! 0
-     X(2) = -X(0)
-     W(0) =  0.555555555555555555555555555555555556_rk ! 5/9
-     W(1) =  0.888888888888888888888888888888888889_rk ! 8/9
-     W(2) =  W(0)
-  case (4) ! p <= 7
-     X(0) = -0.861136311594052575223946488892809506_rk ! sqrt((3+2*sqrt(6/5))/7)
-     X(1) = -0.339981043584856264802665759103244686_rk ! sqrt((3-2*sqrt(6/5))/7)
-     X(2) = -X(1)
-     X(3) = -X(0)
-     W(0) =  0.347854845137453857373063949221999408_rk ! (18-sqrt(30))/36
-     W(1) =  0.652145154862546142626936050778000592_rk ! (18+sqrt(30))/36
-     W(2) =  W(1)
-     W(3) =  W(0)
-  case (5) ! p <= 9
-     X(0) = -0.906179845938663992797626878299392962_rk ! 1/3*sqrt(5+2*sqrt(10/7))
-     X(1) = -0.538469310105683091036314420700208806_rk ! 1/3*sqrt(5-2*sqrt(10/7))
-     X(2) =  0.0_rk                                    ! 0
-     X(3) = -X(1)
-     X(4) = -X(0)
-     W(0) =  0.236926885056189087514264040719917362_rk ! (322-13*sqrt(70))/900
-     W(1) =  0.478628670499366468041291514835638193_rk ! (322+13*sqrt(70))/900
-     W(2) =  0.568888888888888888888888888888888889_rk ! 128/225
-     W(3) =  W(1)
-     W(4) =  W(0)
-  case (6) ! p <= 11
-     X(0) = -0.9324695142031520278123015544939946_rk ! << NumericalDifferentialEquationAnalysis`
-     X(1) = -0.6612093864662645136613995950199053_rk ! GaussianQuadratureWeights(6, -1, 1, 37)
-     X(2) = -0.2386191860831969086305017216807119_rk
-     X(3) = -X(2)
-     X(4) = -X(1)
-     X(5) = -X(0)
-     W(0) =  0.171324492379170345040296142172732894_rk
-     W(1) =  0.360761573048138607569833513837716112_rk
-     W(2) =  0.467913934572691047389870343989550995_rk
-     W(3) =  W(2)
-     W(4) =  W(1)
-     W(5) =  W(0)
-  case (7) ! p <= 13
-     X(0) = -0.9491079123427585245261896840478513_rk ! << NumericalDifferentialEquationAnalysis`
-     X(1) = -0.7415311855993944398638647732807884_rk ! GaussianQuadratureWeights(7, -1, 1, 37)
-     X(2) = -0.4058451513773971669066064120769615_rk
-     X(3) =  0.0_rk
-     X(4) = -X(2)
-     X(5) = -X(1)
-     X(6) = -X(0)
-     W(0) =  0.129484966168869693270611432679082018_rk
-     W(1) =  0.279705391489276667901467771423779582_rk
-     W(2) =  0.381830050505118944950369775488975134_rk
-     W(3) =  0.417959183673469387755102040816326531_rk
-     W(4) =  W(2)
-     W(5) =  W(1)
-     W(6) =  W(0)
-  case (8) ! p <= 15
-     X(0) = -0.9602898564975362316835608685694730_rk ! << NumericalDifferentialEquationAnalysis`
-     X(1) = -0.7966664774136267395915539364758304_rk ! GaussianQuadratureWeights(8, -1, 1, 37)
-     X(2) = -0.5255324099163289858177390491892463_rk
-     X(3) = -0.1834346424956498049394761423601840_rk
-     X(4) = -X(3)
-     X(5) = -X(2)
-     X(6) = -X(1)
-     X(7) = -X(0)
-     W(0) =  0.101228536290376259152531354309962190_rk
-     W(1) =  0.222381034453374470544355994426240884_rk
-     W(2) =  0.313706645877887287337962201986601313_rk
-     W(3) =  0.362683783378361982965150449277195612_rk
-     W(4) =  W(3)
-     W(5) =  W(2)
-     W(6) =  W(1)
-     W(7) =  W(0)
-  case (9) ! p <= 17
-     X(0) = -0.9681602395076260898355762029036729_rk ! << NumericalDifferentialEquationAnalysis`
-     X(1) = -0.8360311073266357942994297880697349_rk ! GaussianQuadratureWeights(9, -1, 1, 37)
-     X(2) = -0.6133714327005903973087020393414742_rk
-     X(3) = -0.3242534234038089290385380146433366_rk
-     X(4) =  0.0_rk
-     X(5) = -X(3)
-     X(6) = -X(2)
-     X(7) = -X(1)
-     X(8) = -X(0)
-     W(0) =  0.081274388361574411971892158110523651_rk
-     W(1) =  0.180648160694857404058472031242912810_rk
-     W(2) =  0.260610696402935462318742869418632850_rk
-     W(3) =  0.312347077040002840068630406584443666_rk
-     W(4) =  0.330239355001259763164525069286974049_rk
-     W(5) =  W(3)
-     W(6) =  W(2)
-     W(7) =  W(1)
-     W(8) =  W(0)
-  case default
-     X = 0.0_rk
-     W = 0.0_rk
-  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)
-  integer, parameter :: rk = IGA_REAL_KIND
-  select case (q)
-  case (2) ! p <= 1
-     X(0) = -1.0_rk
-     X(1) = -X(0)
-     W(0) =  1.0_rk
-     W(1) =  W(0)
-  case (3) ! p <= 3
-     X(0) = -1.0_rk                                    ! -1
-     X(1) =  0.0_rk                                    ! 0
-     X(2) = -X(0)
-     W(0) =  0.333333333333333333333333333333333333_rk ! 1/3
-     W(1) =  1.333333333333333333333333333333333333_rk ! 4/3
-     W(2) =  W(0)
-  case (4) ! p <= 5
-     X(0) = -1.0_rk                                    ! -1
-     X(1) = -0.447213595499957939281834733746255246_rk ! 1/sqrt(5)
-     X(2) = -X(1)
-     X(3) = -X(0)
-     W(0) =  0.166666666666666666666666666666666667_rk ! 1/6
-     W(1) =  0.833333333333333333333333333333333343_rk ! 5/6
-     W(2) =  W(1)
-     W(3) =  W(0)
-  case (5) ! p <= 7
-     X(0) = -1.0_rk                                    ! -1
-     X(1) = -0.654653670707977143798292456246858356_rk ! sqrt(3/7)
-     X(2) =  0.0_rk                                    ! 0
-     X(3) = -X(1)
-     X(4) = -X(0)
-     W(0) =  0.1_rk                                    !  1/10
-     W(1) =  0.544444444444444444444444444444444444_rk ! 49/90
-     W(2) =  0.711111111111111111111111111111111111_rk ! 32/45
-     W(3) =  W(1)
-     W(4) =  W(0)
-  case default
-     X = 0.0_rk
-     W = 0.0_rk
-  end select
-end subroutine IGA_Rule_GaussLobatto