Lisandro Dalcin avatar Lisandro Dalcin committed ca3d144

Better Fortran API for element routines

Comments (0)

Files changed (15)

 #include "petiga.h"
 
+typedef struct { 
+  PetscReal lambda;
+} AppCtx;
+
 EXTERN_C_BEGIN
 extern PetscErrorCode Bratu_Function(IGAPoint,const PetscScalar U[],PetscScalar F[],void *);
 extern PetscErrorCode Bratu_Jacobian(IGAPoint,const PetscScalar U[],PetscScalar J[],void *);
 EXTERN_C_END
 
+EXTERN_C_BEGIN
+extern PetscErrorCode Bratu_IFunction(IGAPoint,PetscReal dt,
+                                      PetscReal a,const PetscScalar *V,
+                                      PetscReal t,const PetscScalar *U,
+                                      PetscScalar *F,void *ctx);
+extern PetscErrorCode Bratu_IJacobian(IGAPoint,PetscReal dt,
+                                      PetscReal a,const PetscScalar *V,
+                                      PetscReal t,const PetscScalar *U,
+                                      PetscScalar *F,void *ctx);
+EXTERN_C_END
 
 #undef __FUNCT__
 #define __FUNCT__ "main"
 
   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+  
+  
+  PetscBool steady = PETSC_TRUE;
+  PetscReal lambda = 6.80;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Bratu Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-lambda","Bratu parameter",__FILE__,lambda,&lambda,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-steady","Steady problem",__FILE__,steady,&steady,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
 
     }
   }
 
-  IGAUserFunction Function = Bratu_Function;
-  IGAUserJacobian Jacobian = Bratu_Jacobian;
-  ierr = IGASetUserFunction(iga,*Function,0);CHKERRQ(ierr);
-  ierr = IGASetUserJacobian(iga,*Jacobian,0);CHKERRQ(ierr);
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
 
   PetscInt dim;
   if (dim < 1) {ierr = IGASetDim(iga,dim=2);CHKERRQ(ierr);}
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
+  AppCtx ctx;
+  ctx.lambda = lambda;
+  if (steady) {
+    ierr = IGASetUserFunction(iga,Bratu_Function,&ctx);CHKERRQ(ierr);
+    ierr = IGASetUserJacobian(iga,Bratu_Jacobian,&ctx);CHKERRQ(ierr);
+  } else {
+    ierr = IGASetUserIFunction(iga,Bratu_IFunction,&ctx);CHKERRQ(ierr);
+    ierr = IGASetUserIJacobian(iga,Bratu_IJacobian,&ctx);CHKERRQ(ierr);
+  }
+
   Vec x;
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
-  SNES snes;
-  ierr = IGACreateSNES(iga,&snes);CHKERRQ(ierr);
-
-  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
-  ierr = SNESSolve(snes,0,x);CHKERRQ(ierr);
-
+  if (steady) {
+    SNES snes;
+    ierr = IGACreateSNES(iga,&snes);CHKERRQ(ierr);
+    ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
+    ierr = SNESSolve(snes,0,x);CHKERRQ(ierr);
+    ierr = SNESDestroy(&snes);CHKERRQ(ierr);
+  } else {
+    TS ts;
+    ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
+    ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr);
+    ierr = TSSetDuration(ts,10000,0.1);CHKERRQ(ierr);
+    ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr);
+    ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+    ierr = TSSolve(ts,x,PETSC_NULL);CHKERRQ(ierr);
+    ierr = TSDestroy(&ts);CHKERRQ(ierr);
+  }
   PetscBool draw = PETSC_FALSE;
   ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
   if (draw && dim < 3) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+  ierr = VecDestroy(&x);CHKERRQ(ierr);
 
-  ierr = VecDestroy(&x);CHKERRQ(ierr);
-  ierr = SNESDestroy(&snes);CHKERRQ(ierr);
   ierr = IGADestroy(&iga);CHKERRQ(ierr);
 
   PetscBool flag = PETSC_FALSE;
 #define scalar real
 #endif
 
-!module Bratu 
-!contains
+module BratuFJ 
 
-integer(kind=IGA_ERRCODE) &
+use PetIGA
+implicit none
+
+type, bind(C) :: IGAUser
+   real(kind=IGA_REAL_KIND) lambda
+end type IGAUser
+
+contains
+
+! --- Steady ---
+
+integer(kind=IGA_ERRCODE_KIND) &
 function Function(p,UU,FF,ctx) result (ierr) &
 bind(C, name="Bratu_Function")
   use PetIGA
   implicit none
-  type(C_PTR), intent(in), value :: ctx
-  type(IGAPoint), intent(in)     :: p
-  scalar (kind=IGA_SCALAR), intent(in)  :: UU(p%nen)
-  scalar (kind=IGA_SCALAR), intent(out) :: FF(p%nen)
-  integer(kind=IGA_INT)    :: a
-  scalar (kind=IGA_SCALAR) :: u, grad_u(p%dim)
-  real   (kind=IGA_REAL)   :: N(p%nen), grad_N(p%dim,p%nen)
-  ierr = IGAPointFormValue(p,UU,u)
-  ierr = IGAPointFormGrad (p,UU,grad_u)
-  ierr = IGAPointFormShapeFuns(p,0,N)
-  ierr = IGAPointFormShapeFuns(p,1,grad_N)
+  type(IGAPoint), intent(in) :: p
+  type(IGAUser),  intent(in) :: ctx
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: FF(p%nen)
+  real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:)
+  scalar (kind=IGA_SCALAR_KIND ) :: u, grad_u(p%dim)
+  integer(kind=IGA_INTEGER_KIND) :: a
+
+  N      => IGA_Shape0(p)
+  grad_N => IGA_Shape1(p)
+
+  u      = IGA_Eval(N,UU)
+  grad_u = IGA_Eval(grad_N,UU)
+
+  u      = IGA_Value(p,UU) ! just for testing
+  grad_u = IGA_Grad(p,UU)  ! just for testing
+
   do a = 1, p%nen
-     FF(a) = dot_product(grad_N(:,a),real(grad_u)) - N(a) * 1 * exp(real(u))
+     FF(a) = dot_product(grad_N(:,a),real(grad_u)) - &
+             N(a) * ctx%lambda * exp(real(u))
   end do
+
   ierr = 0
 end function Function
 
-integer(kind=IGA_ERRCODE) &
+integer(kind=IGA_ERRCODE_KIND) &
 function Jacobian(p,UU,JJ,ctx) result (ierr) &
 bind(C, name="Bratu_Jacobian")
   use PetIGA
   implicit none
-  type(C_PTR), intent(in), value :: ctx
-  type(IGAPoint), intent(in)     :: p
-  scalar (kind=IGA_SCALAR), intent(in)  :: UU(p%nen)
-  scalar (kind=IGA_SCALAR), intent(out) :: JJ(p%nen,p%nen)
-  integer(kind=IGA_INT )    :: a, b
-  scalar (kind=IGA_SCALAR)  :: u
-  real   (kind=IGA_REAL)    :: N(p%nen), grad_N(p%dim,p%nen)
-  ierr = IGAPointFormValue(p,UU,u)
-  ierr = IGAPointFormShapeFuns(p,0,N)
-  ierr = IGAPointFormShapeFuns(p,1,grad_N)
+  type(IGAPoint), intent(in) :: p
+  type(IGAUser),  intent(in) :: ctx
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: JJ(p%nen,p%nen)
+  real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:)
+  scalar (kind=IGA_SCALAR_KIND ) :: u
+  integer(kind=IGA_INTEGER_KIND) :: a, b
+
+  N      => IGA_Shape0(p)
+  grad_N => IGA_Shape1(p)
+
+  u = IGA_Eval(N,UU)
+
+  u = IGA_Value(p,UU)
+
   do a = 1, p%nen
      do b = 1, p%nen
-        JJ(b,a) = dot_product(grad_N(:,a),grad_N(:,b)) - N(a) * N(b) * 1 * exp(real(u))
+        JJ(b,a) = dot_product(grad_N(:,a),grad_N(:,b)) - &
+                  N(a) * N(b) * ctx%lambda  * exp(real(u))
      end do
   end do
+
   ierr = 0
 end function Jacobian
 
-!end module Bratu 
+! --- Time-dependent ---
+
+integer(kind=IGA_ERRCODE_KIND) &
+function IFunction(p,dt,shift,VV,t,UU,FF,ctx) result (ierr) &
+bind(C, name="Bratu_IFunction")
+  use PetIGA
+  implicit none
+  type(IGAPoint), intent(in) :: p
+  type(IGAUser),  intent(in) :: ctx
+  real   (kind=IGA_REAL_KIND   ), intent(in), value :: dt, shift, t
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: VV(p%nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: FF(p%nen)
+
+  real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:)
+  scalar (kind=IGA_SCALAR_KIND ) :: v, u, grad_u(p%dim)
+  integer(kind=IGA_INTEGER_KIND) :: a
+
+  N      => IGA_Shape0(p)
+  grad_N => IGA_Shape1(p)
+
+  v      = IGA_Eval(N,VV)
+  u      = IGA_Eval(N,UU)
+  grad_u = IGA_Eval(grad_N,UU)
+
+  do a = 1, p%nen
+     FF(a) = N(a) * real(v) + &
+             dot_product(grad_N(:,a),real(grad_u)) - &
+             N(a) * ctx%lambda * exp(real(u))
+  end do
+
+  ierr = 0
+end function IFunction
+
+integer(kind=IGA_ERRCODE_KIND) &
+function IJacobian(p,dt,shift,VV,t,UU,JJ,ctx) result (ierr) &
+bind(C, name="Bratu_IJacobian")
+  use PetIGA
+  implicit none
+  type(IGAPoint), intent(in) :: p
+  type(IGAUser),  intent(in) :: ctx
+  real   (kind=IGA_REAL_KIND   ), intent(in), value :: dt, shift, t
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: VV(p%nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)  :: UU(p%nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out) :: JJ(p%nen,p%nen)
+
+  real   (kind=IGA_REAL_KIND   ), pointer :: N(:), grad_N(:,:)
+  scalar (kind=IGA_SCALAR_KIND ) :: v, u, grad_u(p%dim)
+  integer(kind=IGA_INTEGER_KIND) :: a, b
+
+  N      => IGA_Shape0(p)
+  grad_N => IGA_Shape1(p)
+
+  v      = IGA_Eval(N,VV)
+  u      = IGA_Eval(N,UU)
+  grad_u = IGA_Eval(grad_N,UU)
+
+  do a = 1, p%nen
+     do b = 1, p%nen
+        JJ(b,a) = shift * N(a) * N(b) + &
+                  dot_product(grad_N(:,a),grad_N(:,b)) - &
+                  N(a) * N(b) * ctx%lambda  * exp(real(u))
+     end do
+  end do
+
+  ierr = 0
+end function IJacobian
+
+end module BratuFJ
 TARGETS = InputOutput \
-          L2Proj1D L2Proj2D \
-          Poisson1D Poisson2D Poisson3D \
-          L2Projection \
-          AdvectionDiffusion \
-          Laplace \
-          Poisson \
-          Bratu \
-          PatternFormation \
-          CahnHilliard2D \
-          CahnHilliard3D \
-          NavierStokesKorteweg2D \
-          NavierStokesVMS \
-          Elasticity3D \
-          HyperElasticity \
-          PhaseFieldCrystal2D
+	  L2Proj1D L2Proj2D \
+	  Poisson1D Poisson2D Poisson3D \
+	  L2Projection \
+	  AdvectionDiffusion \
+	  Laplace \
+	  Poisson \
+	  Bratu \
+	  PatternFormation \
+	  CahnHilliard2D \
+	  CahnHilliard3D \
+	  NavierStokesKorteweg2D \
+	  NavierStokesVMS \
+	  Elasticity3D \
+	  HyperElasticity \
+	  PhaseFieldCrystal2D
 
 ALL: ${TARGETS}
 clean::
 	${RM} -f $<
 Bratu: Bratu.o BratuFJ.o chkopts
 	${CLINKER} -o $@ $< BratuFJ.o ${PETIGA_LIB}
-	${RM} -f $< BratuFJ.o
+	${RM} -f $< BratuFJ.o bratufj.mod
 AdvectionDiffusion: AdvectionDiffusion.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<
 runex5b_4:
 	-@${MPIEXEC} -n 4 ./PatternFormation ${OPTS} -ts_max_steps 2 -implicit
 runex6a_1:
-	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 1
+	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 1 -lambda 3.0
 runex6a_2:
-	-@${MPIEXEC} -n 2 ./Bratu ${OPTS} -iga_dim 1
+	-@${MPIEXEC} -n 2 ./Bratu ${OPTS} -iga_dim 1 -steady false
 runex6b_1:
 	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2
 runex6b_4:
 	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2
+runex6c_1:
+	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2 -steady false -ts_max_steps 2
+runex6c_4:
+	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady false -ts_max_steps 2
 
 InputOutput := InputOutput.PETSc runex0a_1 runex0a_2 runex0a_3 runex0a_4 runex0a_8 runex0a.rm InputOutput.rm
 L2Proj1D := L2Proj1D.PETSc runex1a_1 runex1a_4 L2Proj1D.rm
 Poisson1D := Poisson1D.PETSc runex2a_1 runex2a_4 Poisson1D.rm
 Poisson2D := Poisson2D.PETSc runex2b_1 runex2b_4 Poisson2D.rm
 Poisson3D := Poisson3D.PETSc runex2c_1 runex2c_4 Poisson3D.rm
-Bratu := Bratu.PETSc runex6a_1 runex6a_2 runex6b_1 runex6b_4 Bratu.rm
+Bratu := Bratu.PETSc runex6a_1 runex6a_2 runex6b_1 runex6b_4 runex6c_1 runex6c_4 Bratu.rm
 CahnHilliard2D := CahnHilliard2D.PETSc runex4_1 runex4_4 CahnHilliard2D.rm
 PatternFormation := PatternFormation.PETSc runex5a_1 runex5a_4 runex5b_1 runex5b_4 PatternFormation.rm
 
 CahnHilliard := $(CahnHilliard2D) $(CahnHilliard3D)
 
 TESTEXAMPLES_C := $(InputOutput) $(L2Proj) $(Poisson) $(Bratu) $(CahnHilliard) $(PatternFormation)
-TESTEXAMPLES_F := 
+TESTEXAMPLES_F :=
 TESTEXAMPLES_FORTRAN:=$(TESTEXAMPLES_F)
 testexamples:
 	-@${OMAKE} tree ACTION=testexamples_C PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR}
 #FFLAGS = ${FFLAGS_STD} -pedantic -Wall -Wextra -fcheck=all
 
 SOURCEH  = ../include/petiga.h petigapc.h petigagrid.h petigapart.h
-SOURCEC  = petiga.c petigareg.c petigaftnc.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 petigaio.c petigagrid.c petigapart.c petigacol.c
+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 petigaio.c petigagrid.c petigapart.c petigacol.c
 SOURCEF1 = petigaftn.F90 petigaval.F90
 SOURCEF2 = petigabsp.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90
 SOURCEF  = ${SOURCEF1} ${SOURCEF2}
-subroutine IGA_Quadrature_1D(&
+pure subroutine IGA_Quadrature_1D(&
      inq,iX,iW,iL,           &
      W,J,X,L)                &
   bind(C, name="IGA_Quadrature_1D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 1
-  integer(kind=IGA_INT ), intent(in),value :: inq
-  real   (kind=IGA_REAL), intent(in)  :: iX(inq), iW(inq), iL
-  real   (kind=IGA_REAL), intent(out) :: W(inq)
-  real   (kind=IGA_REAL), intent(out) :: J(inq)
-  real   (kind=IGA_REAL), intent(out) :: X(dim,inq)
-  real   (kind=IGA_REAL), intent(out) :: L(dim,inq)
-  integer(kind=IGA_INT )  :: iq
+  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), iL
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: L(dim,inq)
+  integer(kind=IGA_INTEGER_KIND)  :: iq
   forall (iq=1:inq)
      !
      W(iq) = iW(iq)
 end subroutine IGA_Quadrature_1D
 
 
-subroutine IGA_BasisFuns_1D(&
+pure subroutine IGA_BasisFuns_1D(&
      order,                 &
      rational,W,            &
      inq,ina,ind,iN,        &
   bind(C, name="IGA_BasisFuns_1D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 1
-  integer(kind=IGA_INT ), intent(in),value :: order
-  integer(kind=IGA_INT ), intent(in),value :: rational
-  integer(kind=IGA_INT ), intent(in),value :: inq, ina, ind
-  real   (kind=IGA_REAL), intent(in)  :: iN(0:ind,ina,inq)
-  real   (kind=IGA_REAL), intent(in)  :: W(dim+1,  ina)
-  real   (kind=IGA_REAL), intent(out) :: N0(       ina,inq)
-  real   (kind=IGA_REAL), intent(out) :: N1(   dim,ina,inq)
-  real   (kind=IGA_REAL), intent(out) :: N2(dim**2,ina,inq)
-  real   (kind=IGA_REAL), intent(out) :: N3(dim**3,ina,inq)
-  integer(kind=IGA_INT ) :: ia,iq
-  integer(kind=IGA_INT ) :: nen
+  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)  :: iN(0:ind,ina,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: W(dim+1,  ina)
+  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(&
 
 contains
 
-subroutine TensorBasisFuns(&
+pure subroutine TensorBasisFuns(&
      ord,&
      ina,ind,iN,&
      N0,N1,N2,N3)
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 1
-  integer(kind=IGA_INT ), intent(in),value :: ord
-  integer(kind=IGA_INT ), intent(in),value :: ina, ind
-  real   (kind=IGA_REAL), intent(in)  :: iN(0:ind,ina)
-  real   (kind=IGA_REAL), intent(out) :: N0(            ina)
-  real   (kind=IGA_REAL), intent(out) :: N1(        dim,ina)
-  real   (kind=IGA_REAL), intent(out) :: N2(    dim,dim,ina)
-  real   (kind=IGA_REAL), intent(out) :: N3(dim,dim,dim,ina)
-  integer(kind=IGA_INT ) :: ia
+  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
   !
   forall (ia=1:ina)
      N0(ia) = iN(0,ia)
 end subroutine IGA_BasisFuns_1D
 
 
-subroutine IGA_ShapeFuns_1D(&
+pure subroutine IGA_ShapeFuns_1D(&
      order,                 &
      nqp,nen,X,             &
      M0,M1,M2,M3,           &
   bind(C, name="IGA_ShapeFuns_1D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 1
-  integer(kind=IGA_INT ), intent(in),value :: order
-  integer(kind=IGA_INT ), intent(in),value :: nqp
-  integer(kind=IGA_INT ), intent(in),value :: nen
-  real   (kind=IGA_REAL), intent(in)    :: X(dim+1,nen)
-  real   (kind=IGA_REAL), intent(in)    :: M0(       nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M1(dim,   nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N0(       nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N1(dim,   nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL), intent(inout) :: DetJac(nqp)
-  real   (kind=IGA_REAL), intent(out)   :: F(dim,dim,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: G(dim,dim,nqp)
-  integer(kind=IGA_INT )  :: q
-  real   (kind=IGA_REAL)  :: DetF
+  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+1,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(inout) :: DetJac(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: F(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G(dim,dim,nqp)
+  integer(kind=IGA_INTEGER_KIND)  :: q
+  real   (kind=IGA_REAL_KIND)     :: DetF
   do q=1,nqp
      call GeometryMap(&
           order,&
-subroutine IGA_Quadrature_2D(&
+pure subroutine IGA_Quadrature_2D(&
      inq,iX,iW,iL,           &
      jnq,jX,jW,jL,           &
      W,J,X,L)                &
   bind(C, name="IGA_Quadrature_2D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 2
-  integer(kind=IGA_INT ), intent(in),value :: inq
-  integer(kind=IGA_INT ), intent(in),value :: jnq
-  real   (kind=IGA_REAL), intent(in)  :: iX(inq), iW(inq), iL
-  real   (kind=IGA_REAL), intent(in)  :: jX(jnq), jW(jnq), jL
-  real   (kind=IGA_REAL), intent(out) :: W(    inq,jnq)
-  real   (kind=IGA_REAL), intent(out) :: J(    inq,jnq)
-  real   (kind=IGA_REAL), intent(out) :: X(dim,inq,jnq)
-  real   (kind=IGA_REAL), intent(out) :: L(dim,inq,jnq)
-  integer(kind=IGA_INT ) :: iq
-  integer(kind=IGA_INT ) :: jq
+  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), iL
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jX(jnq), jW(jnq), jL
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(    inq,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(    inq,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: L(dim,inq,jnq)
+  integer(kind=IGA_INTEGER_KIND)  :: iq
+  integer(kind=IGA_INTEGER_KIND)  :: jq
   forall (iq=1:inq, jq=1:jnq)
      !
      W(iq,jq) = iW(iq) * jW(jq)
 end subroutine IGA_Quadrature_2D
 
 
-subroutine IGA_BasisFuns_2D(&
+pure subroutine IGA_BasisFuns_2D(&
      order,                 &
      rational,W,            &
      inq,ina,ind,iN,        &
   bind(C, name="IGA_BasisFuns_2D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 2
-  integer(kind=IGA_INT ), intent(in),value :: order
-  integer(kind=IGA_INT ), intent(in),value :: rational
-  integer(kind=IGA_INT ), intent(in),value :: inq, ina, ind
-  integer(kind=IGA_INT ), intent(in),value :: jnq, jna, jnd
-  real   (kind=IGA_REAL), intent(in)  :: iN(0:ind,ina,inq)
-  real   (kind=IGA_REAL), intent(in)  :: jN(0:jnd,jna,jnq)
-  real   (kind=IGA_REAL), intent(in)  :: W(dim+1,  ina*jna)
-  real   (kind=IGA_REAL), intent(out) :: N0(       ina*jna,inq,jnq)
-  real   (kind=IGA_REAL), intent(out) :: N1(   dim,ina*jna,inq,jnq)
-  real   (kind=IGA_REAL), intent(out) :: N2(dim**2,ina*jna,inq,jnq)
-  real   (kind=IGA_REAL), intent(out) :: N3(dim**3,ina*jna,inq,jnq)
-  integer(kind=IGA_INT ) :: ia,iq
-  integer(kind=IGA_INT ) :: ja,jq
-  integer(kind=IGA_INT ) :: ka,kq
-  integer(kind=IGA_INT ) :: nen
+  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)  :: iN(0:ind,ina,inq)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jN(0:jnd,jna,jnq)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: W(dim+1,  ina*jna)
+  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
 
 contains
 
-subroutine TensorBasisFuns(&
+pure subroutine TensorBasisFuns(&
      ord,&
      ina,ind,iN,&
      jna,jnd,jN,&
      N0,N1,N2,N3)
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 2
-  integer(kind=IGA_INT ), intent(in),value :: ord
-  integer(kind=IGA_INT ), intent(in),value :: ina, ind
-  integer(kind=IGA_INT ), intent(in),value :: jna, jnd
-  real   (kind=IGA_REAL), intent(in)  :: iN(0:ind,ina)
-  real   (kind=IGA_REAL), intent(in)  :: jN(0:jnd,jna)
-  real   (kind=IGA_REAL), intent(out) :: N0(            ina,jna)
-  real   (kind=IGA_REAL), intent(out) :: N1(        dim,ina,jna)
-  real   (kind=IGA_REAL), intent(out) :: N2(    dim,dim,ina,jna)
-  real   (kind=IGA_REAL), intent(out) :: N3(dim,dim,dim,ina,jna)
-  integer(kind=IGA_INT ) :: ia, ja
+  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
   !
   forall (ia=1:ina, ja=1:jna)
      N0(ia,ja) = iN(0,ia) * jN(0,ja)
 end subroutine IGA_BasisFuns_2D
 
 
-subroutine IGA_ShapeFuns_2D(&
+pure subroutine IGA_ShapeFuns_2D(&
      order,                 &
      nqp,nen,X,             &
      M0,M1,M2,M3,           &
   bind(C, name="IGA_ShapeFuns_2D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 2
-  integer(kind=IGA_INT ), intent(in),value :: order
-  integer(kind=IGA_INT ), intent(in),value :: nqp
-  integer(kind=IGA_INT ), intent(in),value :: nen
-  real   (kind=IGA_REAL), intent(in)    :: X(dim+1,nen)
-  real   (kind=IGA_REAL), intent(in)    :: M0(       nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M1(dim,   nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N0(       nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N1(dim,   nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL), intent(inout) :: DetJac(nqp)
-  real   (kind=IGA_REAL), intent(out)   :: F(dim,dim,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: G(dim,dim,nqp)
-  integer(kind=IGA_INT )  :: q
-  real   (kind=IGA_REAL)  :: DetF
+  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+1,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(inout) :: DetJac(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: F(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G(dim,dim,nqp)
+  integer(kind=IGA_INTEGER_KIND)  :: q
+  real   (kind=IGA_REAL_KIND)     :: DetF
   do q=1,nqp
      call GeometryMap(&
           order,&
-subroutine IGA_Quadrature_3D(&
+pure subroutine IGA_Quadrature_3D(&
      inq,iX,iW,iL,           &
      jnq,jX,jW,jL,           &
      knq,kX,kW,kL,           &
   bind(C, name="IGA_Quadrature_3D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 3
-  integer(kind=IGA_INT ), intent(in),value :: inq
-  integer(kind=IGA_INT ), intent(in),value :: jnq
-  integer(kind=IGA_INT ), intent(in),value :: knq
-  real   (kind=IGA_REAL), intent(in)  :: iX(inq), iW(inq), iL
-  real   (kind=IGA_REAL), intent(in)  :: jX(jnq), jW(jnq), jL
-  real   (kind=IGA_REAL), intent(in)  :: kX(knq), kW(knq), kL
-  real   (kind=IGA_REAL), intent(out) :: W(    inq,jnq,knq)
-  real   (kind=IGA_REAL), intent(out) :: J(    inq,jnq,knq)
-  real   (kind=IGA_REAL), intent(out) :: X(dim,inq,jnq,knq)
-  real   (kind=IGA_REAL), intent(out) :: L(dim,inq,jnq,knq)
-  integer(kind=IGA_INT ) :: iq
-  integer(kind=IGA_INT ) :: jq
-  integer(kind=IGA_INT ) :: kq
+  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), iL
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: jX(jnq), jW(jnq), jL
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: kX(knq), kW(knq), kL
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: W(    inq,jnq,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: J(    inq,jnq,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: X(dim,inq,jnq,knq)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: L(dim,inq,jnq,knq)
+  integer(kind=IGA_INTEGER_KIND) :: iq
+  integer(kind=IGA_INTEGER_KIND) :: jq
+  integer(kind=IGA_INTEGER_KIND) :: kq
   forall (iq=1:inq, jq=1:jnq, kq=1:knq)
      !
      W(iq,jq,kq) = iW(iq) * jW(jq) * kW(kq)
 end subroutine IGA_Quadrature_3D
 
 
-subroutine IGA_BasisFuns_3D(&
+pure subroutine IGA_BasisFuns_3D(&
      order,                 &
      rational,W,            &
      inq,ina,ind,iN,        &
   bind(C, name="IGA_BasisFuns_3D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 3
-  integer(kind=IGA_INT ), intent(in),value :: order
-  integer(kind=IGA_INT ), intent(in),value :: rational
-  integer(kind=IGA_INT ), intent(in),value :: inq, ina, ind
-  integer(kind=IGA_INT ), intent(in),value :: jnq, jna, jnd
-  integer(kind=IGA_INT ), intent(in),value :: knq, kna, knd
-  real   (kind=IGA_REAL), intent(in)  :: iN(0:ind,ina,inq)
-  real   (kind=IGA_REAL), intent(in)  :: jN(0:jnd,jna,jnq)
-  real   (kind=IGA_REAL), intent(in)  :: kN(0:knd,kna,knq)
-  real   (kind=IGA_REAL), intent(in)  :: W(dim+1,  ina*jna*kna)
-  real   (kind=IGA_REAL), intent(out) :: N0(       ina*jna*kna,inq,jnq,knq)
-  real   (kind=IGA_REAL), intent(out) :: N1(   dim,ina*jna*kna,inq,jnq,knq)
-  real   (kind=IGA_REAL), intent(out) :: N2(dim**2,ina*jna*kna,inq,jnq,knq)
-  real   (kind=IGA_REAL), intent(out) :: N3(dim**3,ina*jna*kna,inq,jnq,knq)
-  integer(kind=IGA_INT ) :: ia,iq
-  integer(kind=IGA_INT ) :: ja,jq
-  integer(kind=IGA_INT ) :: ka,kq
-  integer(kind=IGA_INT ) :: nen
+  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)  :: 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(in)  :: W(dim+1,  ina*jna*kna)
+  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
 
 contains
 
-subroutine TensorBasisFuns(&
+pure subroutine TensorBasisFuns(&
      ord,&
      ina,ind,iN,&
      jna,jnd,jN,&
      kna,knd,kN,&
      N0,N1,N2,N3)
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 3
-  integer(kind=IGA_INT ), intent(in),value :: ord
-  integer(kind=IGA_INT ), intent(in),value :: ina, ind
-  integer(kind=IGA_INT ), intent(in),value :: jna, jnd
-  integer(kind=IGA_INT ), intent(in),value :: kna, knd
-  real   (kind=IGA_REAL), intent(in)  :: iN(0:ind,ina)
-  real   (kind=IGA_REAL), intent(in)  :: jN(0:jnd,jna)
-  real   (kind=IGA_REAL), intent(in)  :: kN(0:knd,kna)
-  real   (kind=IGA_REAL), intent(out) :: N0(            ina,jna,kna)
-  real   (kind=IGA_REAL), intent(out) :: N1(        dim,ina,jna,kna)
-  real   (kind=IGA_REAL), intent(out) :: N2(    dim,dim,ina,jna,kna)
-  real   (kind=IGA_REAL), intent(out) :: N3(dim,dim,dim,ina,jna,kna)
-  integer(kind=IGA_INT ) :: ia, ja, ka
+  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
   !
   forall (ia=1:ina, ja=1:jna, ka=1:kna)
      N0(ia,ja,ka) = iN(0,ia) * jN(0,ja) * kN(0,ka)
 end subroutine IGA_BasisFuns_3D
 
 
-subroutine IGA_ShapeFuns_3D(&
+pure subroutine IGA_ShapeFuns_3D(&
      order,                 &
      nqp,nen,X,             &
      M0,M1,M2,M3,           &
   bind(C, name="IGA_ShapeFuns_3D")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), parameter        :: dim = 3
-  integer(kind=IGA_INT ), intent(in),value :: order
-  integer(kind=IGA_INT ), intent(in),value :: nqp
-  integer(kind=IGA_INT ), intent(in),value :: nen
-  real   (kind=IGA_REAL), intent(in)    :: X(dim+1,nen)
-  real   (kind=IGA_REAL), intent(in)    :: M0(       nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M1(dim,   nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL), intent(in)    :: M3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N0(       nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N1(dim,   nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N2(dim**2,nen,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: N3(dim**3,nen,nqp)
-  real   (kind=IGA_REAL), intent(inout) :: DetJac(nqp)
-  real   (kind=IGA_REAL), intent(out)   :: F(dim,dim,nqp)
-  real   (kind=IGA_REAL), intent(out)   :: G(dim,dim,nqp)
-  integer(kind=IGA_INT )  :: q
-  real   (kind=IGA_REAL)  :: DetF
+  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+1,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(inout) :: DetJac(nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: F(dim,dim,nqp)
+  real   (kind=IGA_REAL_KIND   ), intent(out)   :: G(dim,dim,nqp)
+  integer(kind=IGA_INTEGER_KIND)  :: q
+  real   (kind=IGA_REAL_KIND   )  :: DetF
   do q=1,nqp
      call GeometryMap(&
           order,&

src/petigabsp.f90

   interface
      pure subroutine DersBasisFuns(i,uu,p,d,U,ders)
        use PetIGA
-       integer(kind=IGA_INT),  intent(in)  :: i, p, d
-       real   (kind=IGA_REAL), intent(in)  :: uu, U(0:i+p)
-       real   (kind=IGA_REAL), intent(out) :: ders(0:p,0:d)
+       integer(kind=IGA_INTEGER_KIND), intent(in)  :: i, p, d
+       real   (kind=IGA_REAL_KIND   ), intent(in)  :: uu, U(0:i+p)
+       real   (kind=IGA_REAL_KIND   ), intent(out) :: ders(0:p,0:d)
      end subroutine DersBasisFuns
   end interface
-  integer(kind=IGA_INT),  intent(in),value :: i, p, d
-  real   (kind=IGA_REAL), intent(in),value :: uu
-  real   (kind=IGA_REAL), intent(in)       :: U(0:i+p)
-  real   (kind=IGA_REAL), intent(out)      :: N(0:d,0:p)
-  real   (kind=IGA_REAL)  :: ders(0:p,0:d)
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: i, p, d
+  real   (kind=IGA_REAL_KIND   ), intent(in),value :: uu
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: U(0:i+p)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: N(0:d,0:p)
+  real   (kind=IGA_REAL_KIND   )  :: ders(0:p,0:d)
   call DersBasisFuns(i,uu,p,d,U,ders)
   N = transpose(ders)
 end subroutine IGA_DersBasisFuns
 pure subroutine DersBasisFuns(i,uu,p,n,U,ders)
   use PetIGA
   implicit none
-  integer(kind=IGA_INT),  intent(in)  :: i, p, n
-  real   (kind=IGA_REAL), intent(in)  :: uu, U(0:i+p)
-  real   (kind=IGA_REAL), intent(out) :: ders(0:p,0:n)
-  integer(kind=IGA_INT)   :: j, k, r, s1, s2, rk, pk, j1, j2
-  real   (kind=IGA_REAL)  :: saved, temp, d
-  real   (kind=IGA_REAL)  :: left(p), right(p)
-  real   (kind=IGA_REAL)  :: ndu(0:p,0:p), a(0:1,0:p)
+  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.0
   do j = 1, p
      left(j)  = uu - U(i+1-j)

src/petigaftn.F90

 #include "petscconf.h"
 
 #define C_PETSC_ERRORCODE C_INT
+#define C_PETSC_BOOL C_INT
 
 #if defined(PETSC_USE_64BIT_INDICES)
 #define C_PETSC_INT C_LONG_LONG
 #endif
 
 #if defined(PETSC_USE_COMPLEX)
-#define C_PETSC_SCALAR  C_PETSC_COMPLEX
+#define C_PETSC_SCALAR C_PETSC_COMPLEX
 #else
-#define C_PETSC_SCALAR  C_PETSC_REAL
+#define C_PETSC_SCALAR C_PETSC_REAL
 #endif
 
 #if defined(PETSC_USE_COMPLEX)
 #endif
 
 module PetIGA
-
-  use ISO_C_BINDING, only: IGA_ERRCODE => C_PETSC_ERRORCODE
   use ISO_C_BINDING, only: C_PTR
-
-  use ISO_C_BINDING, only: IGA_INT     => C_PETSC_INT
-  use ISO_C_BINDING, only: IGA_REAL    => C_PETSC_REAL
-  use ISO_C_BINDING, only: IGA_SCALAR  => C_PETSC_SCALAR
-  use ISO_C_BINDING, only: IGA_COMPLEX => C_PETSC_COMPLEX
+  use ISO_C_BINDING, only: IGA_ERRCODE_KIND => C_PETSC_ERRORCODE
+  use ISO_C_BINDING, only: IGA_LOGICAL_KIND => C_PETSC_BOOL
+  use ISO_C_BINDING, only: IGA_INTEGER_KIND => C_PETSC_INT
+  use ISO_C_BINDING, only: IGA_REAL_KIND    => C_PETSC_REAL
+  use ISO_C_BINDING, only: IGA_COMPLEX_KIND => C_PETSC_COMPLEX
+  use ISO_C_BINDING, only: IGA_SCALAR_KIND  => C_PETSC_SCALAR
   implicit none
 
   !type, bind(C) :: IGA
   !end type IGAElement
 
   !type, bind(C) :: IGAElement
-  !   integer(kind=IGA_INT) :: private_refct
-  !   integer(kind=IGA_INT) :: private_start(3);
-  !   integer(kind=IGA_INT) :: private_width(3);
-  !   integer(kind=IGA_INT) :: ID(3);
-  !   integer(kind=IGA_INT) :: count
-  !   integer(kind=IGA_INT) :: index
-  !   integer(kind=IGA_INT) :: qnp
-  !   integer(kind=IGA_INT) :: nen
-  !   integer(kind=IGA_INT) :: dof
-  !   integer(kind=IGA_INT) :: dim
-  !   integer(kind=IGA_INT) :: nsd
+  !   integer(kind=IGA_INTEGER_KIND) :: private_refct
+  !   integer(kind=IGA_INTEGER_KIND) :: private_start(3);
+  !   integer(kind=IGA_INTEGER_KIND) :: private_width(3);
+  !   integer(kind=IGA_INTEGER_KIND) :: private_ID(3);
+  !   integer(kind=IGA_INTEGER_KIND) :: count
+  !   integer(kind=IGA_INTEGER_KIND) :: index
+  !   integer(kind=IGA_INTEGER_KIND) :: qnp
+  !   integer(kind=IGA_INTEGER_KIND) :: nen
+  !   integer(kind=IGA_INTEGER_KIND) :: dof
+  !   integer(kind=IGA_INTEGER_KIND) :: dim
+  !   integer(kind=IGA_INTEGER_KIND) :: nsd
   !end type IGAElement
 
   type, bind(C) :: IGAPoint
-     integer(kind=IGA_INT) :: private_refct
-     integer(kind=IGA_INT) :: count
-     integer(kind=IGA_INT) :: index
-     integer(kind=IGA_INT) :: nen
-     integer(kind=IGA_INT) :: dof
-     integer(kind=IGA_INT) :: dim
-     integer(kind=IGA_INT) :: nsd
+     integer(kind=IGA_INTEGER_KIND) :: private_refct
+
+     integer(kind=IGA_INTEGER_KIND) :: count
+     integer(kind=IGA_INTEGER_KIND) :: index
+     integer(kind=IGA_INTEGER_KIND) :: nen
+     integer(kind=IGA_INTEGER_KIND) :: dof
+     integer(kind=IGA_INTEGER_KIND) :: dim
+     integer(kind=IGA_INTEGER_KIND) :: nsd
+
+     type(C_PTR) :: weight
+     type(C_PTR) :: detJac
+     type(C_PTR) :: point
+     type(C_PTR) :: scale
+     type(C_PTR) :: basis(0:3)
+     type(C_PTR) :: gradX(0:1)
+     type(C_PTR) :: shape(0:3)
   end type IGAPoint
 
-  interface
 
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointGetIndex(p,index) bind(C)
-       import; implicit none
-       type(IGAPoint), intent(in) :: p
-       integer(kind=IGA_INT), intent(out) :: index
-     end function IGAPointGetIndex
+  interface IGA_GeomMap
+     module procedure IGA_GeomMap
+  end interface IGA_GeomMap
 
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointGetCount(p,count) bind(C)
-       import; implicit none
-       type(IGAPoint), intent(in) :: p
-       integer(kind=IGA_INT), intent(out) :: count
-     end function IGAPointGetCount
+  interface IGA_GradGeomMap
+     module procedure IGA_GradGeomMap
+  end interface IGA_GradGeomMap
 
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointFormPoint(p,x) bind(C)
-       import; implicit none
-       type(IGAPoint), intent(in) :: p
-       real(kind=IGA_REAL), intent(out) :: x(p%dim)
-     end function IGAPointFormPoint
+  !interface IGA_HessGeomMap
+  !   module procedure IGA_HessGeomMap
+  !end interface IGA_HessGeomMap
 
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointFormGradMap(p,map,inv) bind(C)
-       import; implicit none
-       type(IGAPoint), intent(in) :: p
-       real(kind=IGA_REAL), intent(out) :: map(p%dim,p%dim)
-       real(kind=IGA_REAL), intent(out) :: inv(p%dim,p%dim)
-     end function IGAPointFormGradMap
+  interface IGA_InvGradGeomMap
+     module procedure IGA_InvGradGeomMap
+  end interface IGA_InvGradGeomMap
 
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointFormShapeFuns(p,der,N) bind(C)
-       import; implicit none
-       type(IGAPoint), intent(in) :: p
-       integer(kind=IGA_INT), intent(in),value :: der
-       real   (kind=IGA_REAL), intent(out) :: N(p%dim**der,p%nen)
-     end function IGAPointFormShapeFuns
 
-  end interface
+  interface IGA_Shape0
+     module procedure IGA_Shape0
+  end interface IGA_Shape0
 
-  interface IGAPointFormValue
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointFormValue(p,U,v) bind(C)
-       import; implicit none
-       type(IGAPoint), intent(in) :: p
-       scalar(kind=IGA_SCALAR), intent(in)  :: U(p%dof,p%nen)
-       scalar(kind=IGA_SCALAR), intent(out) :: v(p%dof)
-     end function IGAPointFormValue
-     module procedure IGAPointFormValue_S
-  end interface IGAPointFormValue
+  interface IGA_Shape1
+     module procedure IGA_Shape1
+  end interface IGA_Shape1
 
-  interface IGAPointFormGrad
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointFormGrad(p,U,v) bind(C)
-       import; implicit none
-       type(IGAPoint), intent(in) :: p
-       scalar(kind=IGA_SCALAR), intent(in)  :: U(p%dof,p%nen)
-       scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim,p%dof)
-     end function IGAPointFormGrad
-     module procedure IGAPointFormGrad_S
-  end interface IGAPointFormGrad
+  interface IGA_Shape2
+     module procedure IGA_Shape2
+  end interface IGA_Shape2
 
-  interface IGAPointFormHess
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointFormHess(p,U,v) bind(C)
-       import; implicit none
-       integer(kind=IGA_ERRCODE)     :: ierr
-       type(IGAPoint), intent(in) :: p
-       scalar(kind=IGA_SCALAR), intent(in)  :: U(p%dof,p%nen)
-       scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim,p%dim,p%dof)
-     end function IGAPointFormHess
-     module procedure IGAPointFormHess_S
-  end interface IGAPointFormHess
+  interface IGA_Shape3
+     module procedure IGA_Shape3
+  end interface IGA_Shape3
 
-  interface IGAPointFormDer0
-     module procedure IGAPointFormValue_V
-     module procedure IGAPointFormValue_S
-  end interface IGAPointFormDer0
 
-  interface IGAPointFormDer1
-     module procedure IGAPointFormGrad_V
-     module procedure IGAPointFormGrad_S
-  end interface IGAPointFormDer1
+  interface IGA_Eval
+     module procedure IGA_Shape_Der0_S
+     module procedure IGA_Shape_Der0_V
+     module procedure IGA_Shape_Der1_S
+     module procedure IGA_Shape_Der1_V
+     module procedure IGA_Shape_Der2_S
+     module procedure IGA_Shape_Der2_V
+     module procedure IGA_Shape_Der3_S
+     module procedure IGA_Shape_Der3_V
+  end interface IGA_Eval
 
-  interface IGAPointFormDer2
-     module procedure IGAPointFormHess_V
-     module procedure IGAPointFormHess_S
-  end interface IGAPointFormDer2
+  interface IGA_Value
+     module procedure IGA_Shape_Der0_S
+     module procedure IGA_Shape_Der0_V
+     module procedure IGA_Point_Der0_S
+     module procedure IGA_Point_Der0_V
+  end interface IGA_Value
 
-  interface IGAPointFormDer3
-     integer(kind=IGA_ERRCODE) &
-     function IGAPointFormDer3(p,U,v) bind(C)
-       import; implicit none
-       type(IGAPoint), intent(in) :: p
-       scalar(kind=IGA_SCALAR), intent(in)  :: U(p%dof,p%nen)
-       scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim,p%dim,p%dim,p%dof)
-     end function IGAPointFormDer3
-     module procedure IGAPointFormDer3_S
-  end interface IGAPointFormDer3
+  interface IGA_Grad
+     module procedure IGA_Shape_Der1_S
+     module procedure IGA_Shape_Der1_V
+     module procedure IGA_Point_Der1_S
+     module procedure IGA_Point_Der1_V
+  end interface IGA_Grad
+
+  interface IGA_Hess
+     module procedure IGA_Shape_Der2_S
+     module procedure IGA_Shape_Der2_V
+     module procedure IGA_Point_Der2_S
+     module procedure IGA_Point_Der2_V
+  end interface IGA_Hess
+
+  interface IGA_Der3
+     module procedure IGA_Shape_Der3_S
+     module procedure IGA_Shape_Der3_V
+     module procedure IGA_Point_Der3_S
+     module procedure IGA_Point_Der3_V
+  end interface IGA_Der3
+
+  interface IGA_Div
+     module procedure IGA_Shape_Div
+     module procedure IGA_Point_Div
+  end interface IGA_Div
+
+  interface IGA_Del2
+     module procedure IGA_Shape_Del2_S
+     module procedure IGA_Shape_Del2_V
+     module procedure IGA_Point_Del2_S
+     module procedure IGA_Point_Del2_V
+  end interface IGA_Del2
 
   contains
 
-    integer(kind=IGA_ERRCODE) &
-    function IGAPointFormValue_S(p,U,v) bind(C) result (ierr) 
+    pure function IGA_GeomMap(p) result (X)
       implicit none
       type(IGAPoint), intent(in) :: p
-      scalar(kind=IGA_SCALAR), intent(in)  :: U(p%nen)
-      scalar(kind=IGA_SCALAR), intent(out) :: v
-      scalar(kind=IGA_SCALAR)  :: tmp(1)
-      ierr = IGAPointFormValue(p, reshape(U,(/1,p%nen/)), tmp)
-      v = tmp(1)
-    end function IGAPointFormValue_S
+      real(kind=IGA_REAL_KIND)   :: X(p%nsd)
+      interface
+         pure subroutine IGAPoint_GeomMap(p,X) bind(C)
+           import; implicit none;
+           type(IGAPoint), intent(in) :: p
+           real(kind=IGA_REAL_KIND), intent(out) :: X(p%nsd)
+         end subroutine IGAPoint_GeomMap
+      end interface
+      call IGAPoint_GeomMap(p,X)
+    end function IGA_GeomMap
 
-    integer(kind=IGA_ERRCODE) &
-    function IGAPointFormValue_V(p,U,v) bind(C) result (ierr)
+    pure function IGA_GradGeomMap(p) result (F)
       implicit none
       type(IGAPoint), intent(in) :: p
-      scalar(kind=IGA_SCALAR), intent(in)  :: U(p%dof,p%nen)
-      scalar(kind=IGA_SCALAR), intent(out) :: v(p%dof)
-      ierr = IGAPointFormValue(p,U,v)
-    end function IGAPointFormValue_V
+      real(kind=IGA_REAL_KIND)   :: F(p%dim,p%nsd)
+      interface
+         pure subroutine IGAPoint_GradGeomMap(p,F) bind(C)
+           import; implicit none;
+           type(IGAPoint), intent(in) :: p
+           real(kind=IGA_REAL_KIND), intent(out) :: F(p%dim,p%nsd)
+         end subroutine IGAPoint_GradGeomMap
+      end interface
+      call IGAPoint_GradGeomMap(p,F)
+    end function IGA_GradGeomMap
 
-    integer(kind=IGA_ERRCODE) &
-    function IGAPointFormGrad_S(p,U,v) bind(C) result (ierr)
+    pure function IGA_InvGradGeomMap(p) result (G)
       implicit none
       type(IGAPoint), intent(in) :: p
-      scalar(kind=IGA_SCALAR), intent(in)  :: U(p%nen)
-      scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim)
-      scalar(kind=IGA_SCALAR)  :: tmp(p%dim,1)
-      ierr = IGAPointFormGrad_V(p, reshape(U,(/1,p%nen/)), tmp)
-      v = tmp(:,1)
-    end function IGAPointFormGrad_S
+      real(kind=IGA_REAL_KIND)   :: G(p%nsd,p%dim)
+      interface
+         pure subroutine IGAPoint_InvGradGeomMap(p,G) bind(C)
+           import; implicit none;
+           type(IGAPoint), intent(in) :: p
+           real(kind=IGA_REAL_KIND), intent(out) :: G(p%nsd,p%dim)
+         end subroutine IGAPoint_InvGradGeomMap
+      end interface
+      call IGAPoint_InvGradGeomMap(p,G)
+    end function IGA_InvGradGeomMap
 
-    integer(kind=IGA_ERRCODE) &
-    function IGAPointFormGrad_V(p,U,v) bind(C) result (ierr)
+    pure function IGA_Shape0(p) result(N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
       implicit none
       type(IGAPoint), intent(in) :: p
-      scalar(kind=IGA_SCALAR), intent(in)  :: U(p%dof,p%nen)
-      scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim,p%dof)
-      ierr = IGAPointFormGrad(p,U,v)
-    end function IGAPointFormGrad_V
+      real(kind=IGA_REAL_KIND), pointer :: N(:)
+      call c2f(p%shape(0),N,(/p%nen/))
+    end function IGA_Shape0
 
-    integer(kind=IGA_ERRCODE) &
-    function IGAPointFormHess_S(p,U,v) bind(C) result (ierr)
+    pure function IGA_Shape1(p) result(N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
       implicit none
       type(IGAPoint), intent(in) :: p
-      scalar(kind=IGA_SCALAR), intent(in)  :: U(p%nen)
-      scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim,p%dim)
-      scalar(kind=IGA_SCALAR)  :: tmp(p%dim,p%dim,1)
-      ierr = IGAPointFormHess(p, reshape(U,(/1,p%nen/)), tmp)
-      v = tmp(:,:,1)
-    end function IGAPointFormHess_S
+      real(kind=IGA_REAL_KIND), pointer :: N(:,:)
+      call c2f(p%shape(1),N,(/p%dim,p%nen/))
+    end function IGA_Shape1
 
-    integer(kind=IGA_ERRCODE) &
-    function IGAPointFormHess_V(p,U,v) bind(C) result (ierr)
+    pure function IGA_Shape2(p) result(N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
       implicit none
       type(IGAPoint), intent(in) :: p
-      scalar(kind=IGA_SCALAR), intent(in)  :: U(p%dof,p%nen)
-      scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim,p%dim,p%dof)
-      scalar(kind=IGA_SCALAR)  :: tmp(p%dim,p%dim,1)
-      ierr = IGAPointFormHess(p,U,v)
-    end function IGAPointFormHess_V
+      real(kind=IGA_REAL_KIND), pointer :: N(:,:,:)
+      call c2f(p%shape(2),N,(/p%dim,p%dim,p%nen/))
+    end function IGA_Shape2
 
-    integer(kind=IGA_ERRCODE) &
-    function IGAPointFormDer3_S(p,U,v) bind(C) result (ierr)
+    pure function IGA_Shape3(p) result(N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
       implicit none
       type(IGAPoint), intent(in) :: p
-      scalar(kind=IGA_SCALAR), intent(in)  :: U(p%nen)
-      scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim,p%dim,p%dim)
-      scalar(kind=IGA_SCALAR)  :: tmp(p%dim,p%dim,p%dim,1)
-      ierr = IGAPointFormDer3(p, reshape(U,(/1,p%nen/)), tmp)
-      v = tmp(:,:,:,1)
-    end function IGAPointFormDer3_S
+      real(kind=IGA_REAL_KIND), pointer :: N(:,:,:,:)
+      call c2f(p%shape(3),N,(/p%dim,p%dim,p%dim,p%nen/))
+    end function IGA_Shape3
 
-    integer(kind=IGA_ERRCODE) &
-    function IGAPointFormDer3_V(p,U,v) bind(C) result (ierr)
+#define DIM size(N,1)
+#define DOF size(U,1)
+
+    pure function IGA_Shape_Der0_S(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:) ! nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:) ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V
+      ! V = dot_product(N,U)
+      integer a
+      V = 0
+      do a = 1, size(U,1) ! nen
+         V = V + N(a) * U(a)
+      end do
+    end function IGA_Shape_Der0_S
+
+    pure function IGA_Shape_Der0_V(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:)   ! nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DOF)            ! dof
+      ! V = matmul(N,transpose(U))
+      integer a
+      V = 0
+      do a = 1, size(U,2) ! nen
+         V = V + N(a) * U(:,a)
+      end do
+    end function IGA_Shape_Der0_V
+
+    pure function IGA_Shape_Der1_S(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:) ! dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:)   ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM)            ! dim
+      !V = matmul(N,U)
+      integer a
+      V = 0
+      do a = 1, size(U,1) ! nen
+         V(:) = V(:) + N(:,a) * U(a)
+      end do
+    end function IGA_Shape_Der1_S
+
+    pure function IGA_Shape_Der1_V(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:) ! dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DOF)        ! dim,dof
+      ! V = matmul(N,transpose(U))
+      integer a, c
+      V = 0
+      do a = 1, size(U,2) ! nen
+         do c = 1, size(U,1) ! dof
+            V(:,c) = V(:,c) + N(:,a) * U(c,a)
+         end do
+      end do
+    end function IGA_Shape_Der1_V
+
+    pure function IGA_Shape_Der2_S(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:,:) ! dim,dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:)     ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM)          ! dim,dim
+      integer a
+      V = 0
+      do a = 1, size(U,1) ! nen
+         V(:,:) = V(:,:) + N(:,:,a) * U(a)
+      end do
+    end function IGA_Shape_Der2_S
+
+    pure function IGA_Shape_Der2_V(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:,:) ! dim,dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:)   ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DOF)      ! dim,dim,dof
+      integer a, c
+      V = 0
+      do a = 1, size(U,2) ! nen
+         do c = 1, size(U,1) ! dof
+            V(:,:,c) = V(:,:,c) + N(:,:,a) * U(c,a)
+         end do
+      end do
+    end function IGA_Shape_Der2_V
+
+    pure function IGA_Shape_Der3_S(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:,:,:) ! dim,dim,dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:)       ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DIM)        ! dim,dim,dim
+      integer a
+      V = 0
+      do a = 1, size(U,1) ! nen
+         V(:,:,:) = V(:,:,:) + N(:,:,:,a) * U(a)
+      end do
+    end function IGA_Shape_Der3_S
+
+    pure function IGA_Shape_Der3_V(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:,:,:) ! dim,dim,dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:)     ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DIM,DOF)    ! dim,dim,dim,dof
+      integer a, c
+      V = 0
+      do a = 1, size(U,2) ! nen
+         do c = 1, size(U,1) ! dof
+            V(:,:,:,c) = V(:,:,:,c) + N(:,:,:,a) * U(c,a)
+         end do
+      end do
+    end function IGA_Shape_Der3_V
+
+    pure function IGA_Shape_Div(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:) ! dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dim,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V
+      integer a, c, i
+      V = 0
+      do a = 1, size(U,2) ! nen
+         do i = 1, size(N,1) ! dim
+            V = V + N(i,a) * U(i,a)
+         end do
+      end do
+    end function IGA_Shape_Div
+
+    pure function IGA_Shape_Del2_S(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:,:) ! dim,dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:)     ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V
+      integer a, c, i
+      V = 0
+      do a = 1, size(U,1) ! nen
+         do i = 1, size(N,1) ! dim
+            V = V + N(i,i,a) * U(a)
+         end do
+      end do
+    end function IGA_Shape_Del2_S
+
+    pure function IGA_Shape_Del2_V(N,U) result (V)
+      implicit none
+      real   (kind=IGA_REAL_KIND  ), intent(in) :: N(:,:,:) ! dim,dim,nen
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:)   ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DOF)
+      integer a, c, i
+      V = 0
+      do a = 1, size(U,2) ! nen
+         do c = 1, size(U,1) ! dof
+            do i = 1, size(N,1) ! dim
+               V(c) = V(c) + N(i,i,a) * U(c,a)
+            end do
+         end do
+      end do
+    end function IGA_Shape_Del2_V
+
+#undef DIM
+#undef DOF
+
+#define DIM p%dim
+#define DOF size(U,1)
+
+    pure function IGA_Point_Der0_S(p,U) result(V)
       implicit none
       type(IGAPoint), intent(in) :: p
-      scalar(kind=IGA_SCALAR), intent(in)  :: U(p%dof,p%nen)
-      scalar(kind=IGA_SCALAR), intent(out) :: v(p%dim,p%dim,p%dim,p%dof)
-      ierr = IGAPointFormDer3(p,U,v)
-    end function IGAPointFormDer3_V
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:) ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V               !
+      V = IGA_Eval(IGA_Shape0(p),U)
+    end function IGA_Point_Der0_S
+
+    pure function IGA_Point_Der0_V(p,U) result(V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DOF)            ! dim,dof
+      V = IGA_Eval(IGA_Shape0(p),U)
+    end function IGA_Point_Der0_V
+
+    pure function IGA_Point_Der1_S(p,U) result(V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:) ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM)          ! dim
+      V = IGA_Eval(IGA_Shape1(p),U)
+    end function IGA_Point_Der1_S
+
+    pure function IGA_Point_Der1_V(p,U) result(V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DOF)        ! dim,dof
+      V = IGA_Eval(IGA_Shape1(p),U)
+    end function IGA_Point_Der1_V
+
+    pure function IGA_Point_Der2_S(p,U) result(V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:) ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM)      ! dim,dim
+      V = IGA_Eval(IGA_Shape2(p),U)
+    end function IGA_Point_Der2_S
+
+    pure function IGA_Point_Der2_V(p,U) result(V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DOF)    ! dim,dim,dof
+      V = IGA_Eval(IGA_Shape2(p),U)
+    end function IGA_Point_Der2_V
+
+    pure function IGA_Point_Der3_S(p,U) result(V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:) ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DIM)  ! dim,dim,dim
+      V = IGA_Eval(IGA_Shape3(p),U)
+    end function IGA_Point_Der3_S
+
+    pure function IGA_Point_Der3_V(p,U) result(V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:)  ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DIM,DOF) ! dim,dim,dof
+      V = IGA_Eval(IGA_Shape3(p),U)
+    end function IGA_Point_Der3_V
+
+    pure function IGA_Point_Div(p,U) result (V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dim,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DIM)            ! dim
+      V = IGA_Shape_Div(IGA_Shape1(p),U)
+    end function IGA_Point_Div
+
+    pure function IGA_Point_Del2_S(p,U) result (V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:) ! nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V
+      V = IGA_Shape_Del2_S(IGA_Shape2(p),U)
+    end function IGA_Point_Del2_S
+
+    pure function IGA_Point_Del2_V(p,U) result (V)
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dof,nen
+      scalar (kind=IGA_SCALAR_KIND)  :: V(DOF)
+      V = IGA_Shape_Del2_V(IGA_Shape2(p),U)
+    end function IGA_Point_Del2_V
+
+#undef DIM
+#undef DOF
 
 end module PetIGA

src/petigaftnc.c

-#include "petiga.h"
-
-EXTERN_C_BEGIN
-
-PetscErrorCode igapointgetindex (IGAPoint p,PetscInt *i) {return IGAPointGetIndex(p,i);}
-PetscErrorCode igapointgetcount (IGAPoint p,PetscInt *c) {return IGAPointGetCount(p,c);}
-
-PetscErrorCode igapointformpoint     (IGAPoint p,PetscReal x[])               {return IGAPointFormPoint     (p,x);  }
-PetscErrorCode igapointformgradmap   (IGAPoint p,PetscReal m[],PetscReal i[]) {return IGAPointFormGradMap   (p,m,i);}
-PetscErrorCode igapointformshapefuns (IGAPoint p,PetscInt d,PetscReal N[])    {return IGAPointFormShapeFuns (p,d,N);}
-
-PetscErrorCode igapointformvalue (IGAPoint p,const PetscScalar U[],PetscScalar u[]) {return IGAPointFormValue(p,U,u);}
-PetscErrorCode igapointformgrad  (IGAPoint p,const PetscScalar U[],PetscScalar u[]) {return IGAPointFormGrad (p,U,u);}
-PetscErrorCode igapointformhess  (IGAPoint p,const PetscScalar U[],PetscScalar u[]) {return IGAPointFormHess (p,U,u);}
-PetscErrorCode igapointformdel2  (IGAPoint p,const PetscScalar U[],PetscScalar u[]) {return IGAPointFormDel2 (p,U,u);}
-PetscErrorCode igapointformder3  (IGAPoint p,const PetscScalar U[],PetscScalar u[]) {return IGAPointFormDer3 (p,U,u);}
-
-EXTERN_C_END

src/petigageo.f90.in

      DetG,Grad,InvG)
   use PetIGA
   implicit none
-  !integer(kind=IGA_INT), parameter   :: dim = 1,2,3
-  integer(kind=IGA_INT ), intent(in)  :: order
-  integer(kind=IGA_INT ), intent(in)  :: nen
-  real   (kind=IGA_REAL), intent(in)  :: X(dim,nen)
-  real   (kind=IGA_REAL), intent(in)  :: N0(            nen)
-  real   (kind=IGA_REAL), intent(in)  :: N1(        dim,nen)
-  real   (kind=IGA_REAL), intent(in)  :: N2(    dim,dim,nen)
-  real   (kind=IGA_REAL), intent(in)  :: N3(dim,dim,dim,nen)
-  real   (kind=IGA_REAL), intent(out) :: R0(            nen)
-  real   (kind=IGA_REAL), intent(out) :: R1(        dim,nen)
-  real   (kind=IGA_REAL), intent(out) :: R2(    dim,dim,nen)
-  real   (kind=IGA_REAL), intent(out) :: R3(dim,dim,dim,nen)
-  real   (kind=IGA_REAL), intent(out) :: DetG
-  real   (kind=IGA_REAL), intent(out) :: Grad(dim,dim)
-  real   (kind=IGA_REAL), intent(out) :: InvG(dim,dim)
+  !integer(kind=IGA_INTEGER_KIND),parameter   :: dim = 1,2,3
+  integer(kind=IGA_INTEGER_KIND), intent(in)  :: order
+  integer(kind=IGA_INTEGER_KIND), intent(in)  :: nen
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: X(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: N0(            nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: N1(        dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: N2(    dim,dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)  :: N3(dim,dim,dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: R0(            nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: R1(        dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: R2(    dim,dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: R3(dim,dim,dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: DetG
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: Grad(dim,dim)
+  real   (kind=IGA_REAL_KIND   ), intent(out) :: InvG(dim,dim)
 
-  integer(kind=IGA_INT ) :: idx
-  integer(kind=IGA_INT ) :: i, j, k, l
-  integer(kind=IGA_INT ) :: a, b, c, d
-  real   (kind=IGA_REAL) :: X1(dim,dim)
-  real   (kind=IGA_REAL) :: X2(dim,dim,dim)
-  real   (kind=IGA_REAL) :: X3(dim,dim,dim,dim)
-  real   (kind=IGA_REAL) :: E1(dim,dim)
-  real   (kind=IGA_REAL) :: E2(dim,dim,dim)
-  real   (kind=IGA_REAL) :: E3(dim,dim,dim,dim)
+  integer(kind=IGA_INTEGER_KIND)  :: idx
+  integer(kind=IGA_INTEGER_KIND)  :: i, j, k, l
+  integer(kind=IGA_INTEGER_KIND)  :: a, b, c, d
+  real   (kind=IGA_REAL_KIND   )  :: X1(dim,dim)
+  real   (kind=IGA_REAL_KIND   )  :: X2(dim,dim,dim)
+  real   (kind=IGA_REAL_KIND   )  :: X3(dim,dim,dim,dim)
+  real   (kind=IGA_REAL_KIND   )  :: E1(dim,dim)
+  real   (kind=IGA_REAL_KIND   )  :: E2(dim,dim,dim)
+  real   (kind=IGA_REAL_KIND   )  :: E3(dim,dim,dim,dim)
 
   ! gradient of the mapping
   Grad = matmul(N1,transpose(X))

src/petigainv.f90.in

 
 pure function Determinant(dim,A) result (detA)
   use PetIGA; implicit none
-  integer(kind=IGA_INT ), intent(in) :: dim
-  real   (kind=IGA_REAL), intent(in) :: A(dim,dim)
-  real   (kind=IGA_REAL)  :: detA
+  integer(kind=IGA_INTEGER_KIND), intent(in) :: dim
+  real   (kind=IGA_REAL_KIND   ), intent(in) :: A(dim,dim)
+  real   (kind=IGA_REAL_KIND   )  :: detA
   select case (dim)
   case (1)
      detA = A(1,1)
 
 pure function Inverse(dim,detA,A) result (invA)
   use PetIGA; implicit none
-  integer(kind=IGA_INT ), intent(in) :: dim
-  real   (kind=IGA_REAL), intent(in) :: detA
-  real   (kind=IGA_REAL), intent(in) :: A(dim,dim)
-  real   (kind=IGA_REAL)  :: invA(dim,dim)
+  integer(kind=IGA_INTEGER_KIND), intent(in) :: dim
+  real   (kind=IGA_REAL_KIND   ), intent(in) :: detA
+  real   (kind=IGA_REAL_KIND   ), intent(in) :: A(dim,dim)
+  real   (kind=IGA_REAL_KIND   )  :: invA(dim,dim)
   select case (dim)
   case (1)
      invA = 1.0/detA

src/petigapoint.c

 }
 
 EXTERN_C_BEGIN
-extern void IGA_GetGeomMap (PetscInt nen,PetscInt nsd,
-                            const PetscReal N[],const PetscReal C[],PetscReal X[]);
-extern void IGA_GetGradMap (PetscInt nen,PetscInt nsd,PetscInt dim,
-                            const PetscReal N[],const PetscReal C[],PetscReal F[]);
-extern void IGA_GetGradMapI(PetscInt nen,PetscInt nsd,PetscInt dim,
-                            const PetscReal N[],const PetscReal C[],PetscReal G[]);
+extern void IGA_GetGeomMap     (PetscInt nen,PetscInt nsd,
+                                const PetscReal N[],const PetscReal C[],PetscReal X[]);
+extern void IGA_GetGradGeomMap (PetscInt nen,PetscInt nsd,PetscInt dim,
+                                const PetscReal N[],const PetscReal C[],PetscReal F[]);
+extern void IGA_GetGradGeomMapI(PetscInt nen,PetscInt nsd,PetscInt dim,
+                                const PetscReal N[],const PetscReal C[],PetscReal G[]);
 EXTERN_C_END
 
+#undef  __FUNCT__
+#define __FUNCT__ "IGAPointFormGeomMap"
+PetscErrorCode IGAPointFormGeomMap(IGAPoint p,PetscReal x[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(p,1);
+  PetscValidRealPointer(x,2);
+  if (p->parent->geometry) {
+    const PetscReal *X = p->parent->geometryX;
+    IGA_GetGeomMap(p->nen,p->nsd,p->shape[0],X,x);
+  } else {
+    PetscInt i,dim = p->dim;
+    const PetscReal *X = p->point;
+    for (i=0; i<dim; i++) x[i] = X[i];
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAPointFormGradGeomMap"
+PetscErrorCode IGAPointFormGradGeomMap(IGAPoint p,PetscReal F[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(p,1);
+  PetscValidRealPointer(F,2);
+  if (p->parent->geometry) {
+    PetscInt a,dim = p->dim;
+    PetscInt i,nsd = p->nsd;
+    const PetscReal *L = p->scale;
+    if (dim == nsd) {
+      (void)PetscMemcpy(F,p->gradX[0],nsd*dim*sizeof(PetscReal));
+    } else {
+      const PetscReal *X = p->parent->geometryX;
+      IGA_GetGradGeomMap(p->nen,nsd,dim,p->basis[1],X,F);
+    }
+    for (i=0; i<nsd; i++)
+      for (a=0; a<dim; a++)
+        F[i*dim+a] *= L[a];
+  } else {
+    PetscInt i,dim = p->dim;
+    const PetscReal *L = p->scale;
+    (void)PetscMemzero(F,dim*dim*sizeof(PetscReal));
+    for (i=0; i<dim; i++) F[i*(dim+1)] = L[i];
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAPointFormInvGradGeomMap"
+PetscErrorCode IGAPointFormInvGradGeomMap(IGAPoint p,PetscReal G[])
+{
+  PetscFunctionBegin;
+  PetscValidPointer(p,1);
+  PetscValidRealPointer(G,2);
+  if (p->parent->geometry) {
+    PetscInt a,dim = p->dim;
+    PetscInt i,nsd = p->nsd;
+    const PetscReal *L = p->scale;
+    if (dim == nsd) {
+      (void)PetscMemcpy(G,p->gradX[1],dim*nsd*sizeof(PetscReal));
+    } else {
+      const PetscReal *X = p->parent->geometryX;
+      IGA_GetGradGeomMapI(p->nen,nsd,dim,p->basis[1],X,G);
+    }
+    for (a=0; a<dim; a++)
+      for (i=0; i<nsd; i++)
+        G[a*dim+i] /= L[a];
+  } else {
+    PetscInt i,dim = p->dim;
+    const PetscReal *L = p->scale;
+    (void)PetscMemzero(G,dim*dim*sizeof(PetscReal));
+    for (i=0; i<dim; i++) G[i*(dim+1)] = 1/L[i];
+  }
+  PetscFunctionReturn(0);
+}
+
+EXTERN_C_BEGIN
+void igapoint_geommap (IGAPoint p,PetscReal x[]) {(void)IGAPointFormGeomMap(p,x);}
+void igapoint_gradgeommap (IGAPoint p,PetscReal f[]) {(void)IGAPointFormGradGeomMap(p,f);}
+void igapoint_invgradgeommap(IGAPoint p,PetscReal g[]) {(void)IGAPointFormInvGradGeomMap(p,g);}
+EXTERN_C_END
+
+
 EXTERN_C_BEGIN
 extern void IGA_GetValue(PetscInt nen,PetscInt dof,const PetscReal N[],
                          const PetscScalar U[],PetscScalar u[]);
 #undef  __FUNCT__
 #define __FUNCT__ "IGAPointFormPoint"
 PetscErrorCode IGAPointFormPoint(IGAPoint p,PetscReal x[])
-{
-  PetscFunctionBegin;
-  PetscValidPointer(p,1);
-  PetscValidRealPointer(x,2);
-  if (p->parent->geometry) {
-    const PetscReal *X = p->parent->geometryX;
-    IGA_GetGeomMap(p->nen,p->nsd,p->shape[0],X,x);
-  } else {
-    PetscInt i,dim = p->dim;
-    const PetscReal *X = p->point;
-    for (i=0; i<dim; i++) x[i] = X[i];
-  }
-  PetscFunctionReturn(0);
-}
+{ return IGAPointFormGeomMap(p,x); }
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAPointFormGradMap"
 PetscErrorCode IGAPointFormGradMap(IGAPoint p,PetscReal map[],PetscReal inv[])
 {
-  PetscReal *F = map;
-  PetscReal *G = inv;
+  PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(p,1);
-  if(F) PetscValidRealPointer(F,2);
-  if(G) PetscValidRealPointer(G,3);
-  if (p->parent->geometry) {
-    PetscInt a,dim = p->dim;
-    PetscInt i,nsd = p->nsd;
-    const PetscReal *L = p->scale;
-    if (dim == nsd) {
-      if (F) {(void)PetscMemcpy(F,p->gradX[0],nsd*dim*sizeof(PetscReal));}
-      if (G) {(void)PetscMemcpy(G,p->gradX[1],dim*nsd*sizeof(PetscReal));}
-    } else {
-      const PetscReal *X = p->parent->geometryX;
-      if (F) {IGA_GetGradMap (p->nen,nsd,dim,p->basis[1],X,F);}
-      if (G) {IGA_GetGradMapI(p->nen,nsd,dim,p->basis[1],X,G);}
-    }
-    if (F)
-      for (i=0; i<nsd; i++)
-        for (a=0; a<dim; a++)
-          F[i*dim+a] *= L[a];
-    if (G)
-      for (a=0; a<dim; a++)
-        for (i=0; i<nsd; i++)
-          G[a*dim+i] /= L[a];
-  } else {
-    PetscInt i,dim = p->dim;
-    const PetscReal *L = p->scale;
-    if (F) {
-      (void)PetscMemzero(F,sizeof(PetscReal)*dim*dim);
-      for (i=0; i<dim; i++) F[i*(dim+1)] = L[i];
-    }
-    if (G) {
-      (void)PetscMemzero(G,sizeof(PetscReal)*dim*dim);
-      for (i=0; i<dim; i++) G[i*(dim+1)] = 1/L[i];
-    }
-  }
+  if(map) {ierr = IGAPointFormGradGeomMap(p,map);CHKERRQ(ierr);}
+  if(inv) {ierr = IGAPointFormInvGradGeomMap(p,inv);CHKERRQ(ierr);}
   PetscFunctionReturn(0);
 }
 

src/petigarat.f90.in

 
 pure subroutine Rationalize(ord,nen,W,R0,R1,R2,R3)
   implicit none
-  !integer(kind=IGA_INT), parameter     :: dim = 1,2,3
-  integer(kind=IGA_INT ), intent(in)    :: ord, nen
-  real   (kind=IGA_REAL), intent(in)    :: W(nen)
-  real   (kind=IGA_REAL), intent(inout) :: R0(            nen)
-  real   (kind=IGA_REAL), intent(inout) :: R1(        dim,nen)
-  real   (kind=IGA_REAL), intent(inout) :: R2(    dim,dim,nen)
-  real   (kind=IGA_REAL), intent(inout) :: R3(dim,dim,dim,nen)
+  !integer(kind=IGA_INTEGER_KIND),parameter     :: dim = 1,2,3
+  integer(kind=IGA_INTEGER_KIND), intent(in)    :: ord, 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)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R2(    dim,dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(inout) :: R3(dim,dim,dim,nen)
   !
-  integer(kind=IGA_INT ) :: a, i, j, k
-  real   (kind=IGA_REAL) :: W0, W1(dim), W2(dim,dim), W3(dim,dim,dim)
+  integer(kind=IGA_INTEGER_KIND)  :: a, i, j, k
+  real   (kind=IGA_REAL_KIND   )  :: W0, W1(dim), W2(dim,dim), W3(dim,dim,dim)
   !
   forall (a=1:nen)
      R0(a) = W(a) * R0(a)

src/petigaval.F90

   bind(C, name="IGA_GetGeomMap")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), intent(in),value :: nen,nsd
-  real   (kind=IGA_REAL), intent(in)       :: N(nen)
-  real   (kind=IGA_REAL), intent(in)       :: C(nsd,nen)
-  real   (kind=IGA_REAL), intent(out)      :: X(nsd)
-  !integer(kind=IGA_INT )  :: a
-  !X = 0
-  !do a = 1, nen
-  !   X = X + N(a) * C(:,a)
-  !end do
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,nsd
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: C(nsd,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: X(nsd)
   X = matmul(C,N)
 end subroutine IGA_GetGeomMap
 
-subroutine IGA_GetGradMap(nen,nsd,dim,N,C,F) &
-  bind(C, name="IGA_GetGradMap")
+subroutine IGA_GetGradGeomMap(nen,nsd,dim,N,C,F) &
+  bind(C, name="IGA_GetGradGeomMap")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), intent(in),value :: nen,nsd,dim
-  real   (kind=IGA_REAL), intent(in)       :: N(dim,nen)
-  real   (kind=IGA_REAL), intent(in)       :: C(nsd,nen)
-  real   (kind=IGA_REAL), intent(out)      :: F(dim,nsd)
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,nsd,dim
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: C(nsd,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: F(dim,nsd)
   F = matmul(N,transpose(C))
-end subroutine IGA_GetGradMap
+end subroutine IGA_GetGradGeomMap
 
-subroutine IGA_GetGradMapI(nen,nsd,dim,N,C,G) &
-  bind(C, name="IGA_GetGradMapI")
+subroutine IGA_GetGradGeomMapI(nen,nsd,dim,N,C,G) &
+  bind(C, name="IGA_GetGradGeomMapI")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT ), intent(in),value :: nen,nsd,dim
-  real   (kind=IGA_REAL), intent(in)       :: N(dim,nen)
-  real   (kind=IGA_REAL), intent(in)       :: C(nsd,nen)
-  real   (kind=IGA_REAL), intent(out)      :: G(nsd,dim)
-  real   (kind=IGA_REAL)  :: F(dim,nsd)
-  real   (kind=IGA_REAL)  :: M(nsd,nsd), invM(nsd,nsd)
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,nsd,dim
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: C(nsd,nen)
+  real   (kind=IGA_REAL_KIND   ), intent(out)      :: G(nsd,dim)
+  real   (kind=IGA_REAL_KIND   )  :: F(dim,nsd)
+  real   (kind=IGA_REAL_KIND   )  :: M(nsd,nsd), invM(nsd,nsd)
   F = matmul(N,transpose(C))
   M = matmul(transpose(F),F)
   invM = Inverse(nsd,Determinant(nsd,M),M)
   G = matmul(invM,transpose(F))
 contains
 include 'petigainv.f90.in'
-end subroutine IGA_GetGradMapI
+end subroutine IGA_GetGradGeomMapI
 
 
 subroutine IGA_GetValue(nen,dof,N,U,V) &
   bind(C, name="IGA_GetValue")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof
-  real   (kind=IGA_REAL  ), intent(in)       :: N(nen)
-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dof)
-  integer(kind=IGA_INT   )  :: a, i
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)       :: U(dof,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out)      :: V(dof)
+  integer(kind=IGA_INTEGER_KIND)  :: a, i
   ! V = matmul(N,transpose(U))
   V = 0
   do a = 1, nen
   bind(C, name="IGA_GetGrad")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim
-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim,nen)
-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim,dof)
-  integer(kind=IGA_INT   )  :: a, c
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof,dim
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)       :: U(dof,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out)      :: V(dim,dof)
+  integer(kind=IGA_INTEGER_KIND)  :: a, c
   ! V = matmul(N,transpose(U))
   V = 0
   do a = 1, nen
   bind(C, name="IGA_GetHess")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim
-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim*dim,nen)
-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim*dim,dof)
-  integer(kind=IGA_INT   )  :: a, i
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof,dim
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim*dim,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)       :: U(dof,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out)      :: V(dim*dim,dof)
+  integer(kind=IGA_INTEGER_KIND)  :: a, i
   ! V = matmul(N,transpose(U))
   V = 0
   do a = 1, nen
   bind(C, name="IGA_GetDel2")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim
-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim,dim,nen)
-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dof)
-  integer(kind=IGA_INT   )  :: a, c, i
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof,dim
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim,dim,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)       :: U(dof,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out)      :: V(dof)
+  integer(kind=IGA_INTEGER_KIND)  :: a, c, i
   V = 0
   do a = 1, nen
      do c = 1, dof
      bind(C, name="IGA_GetDer3")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof,dim
-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim*dim*dim,nen)
-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim*dim*dim,dof)
-  integer(kind=IGA_INT   )  :: a, i
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof,dim
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim*dim*dim,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)       :: U(dof,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out)      :: V(dim*dim*dim,dof)
+  integer(kind=IGA_INTEGER_KIND)  :: a, i
   ! V = matmul(N,transpose(U))
   V = 0
   do a = 1, nen
 !  bind(C, name="IGA_GetDerivative")
 !  use PetIGA
 !  implicit none
-!  integer(kind=IGA_INT   ), intent(in),value :: nen,dof
-!  integer(kind=IGA_INT   ), intent(in),value :: dim,der
-!  real   (kind=IGA_REAL  ), intent(in)       :: N(dim**der,nen)
-!  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-!  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim**der,dof)
-!  integer(kind=IGA_INT   )  :: a, i
+!  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof
+!  integer(kind=IGA_INTEGER_KIND), intent(in),value :: dim,der
+!  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim**der,nen)
+!  scalar (kind=IGA_SCALAR_KIND ), intent(in)       :: U(dof,nen)
+!  scalar (kind=IGA_SCALAR_KIND ), intent(out)      :: V(dim**der,dof)
+!  integer(kind=IGA_INTEGER_KIND)  :: a, i
 !  ! V = matmul(N,transpose(U))
 !  V = 0
 !  do a = 1, nen
   bind(C, name="IGA_Interpolate")
   use PetIGA
   implicit none
-  integer(kind=IGA_INT   ), intent(in),value :: nen,dof
-  integer(kind=IGA_INT   ), intent(in),value :: dim,der
-  real   (kind=IGA_REAL  ), intent(in)       :: N(dim**der,nen)
-  scalar (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  scalar (kind=IGA_SCALAR), intent(out)      :: V(dim**der,dof)
-  integer(kind=IGA_INT   )  :: a, i
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: nen,dof
+  integer(kind=IGA_INTEGER_KIND), intent(in),value :: dim,der
+  real   (kind=IGA_REAL_KIND   ), intent(in)       :: N(dim**der,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(in)       :: U(dof,nen)
+  scalar (kind=IGA_SCALAR_KIND ), intent(out)      :: V(dim**der,dof)
+  integer(kind=IGA_INTEGER_KIND)  :: a, i
   ! V = matmul(N,transpose(U))
   V = 0
   do a = 1, nen
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.