Commits

Lisandro Dalcin committed 5ad48db

First steps towards Fortran support for element routines

  • Participants
  • Parent commits 9bb9f2d

Comments (0)

Files changed (9)

File demo/Bratu.c

+#include "petiga.h"
+
+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
+
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  IGABoundary bnd;
+  PetscInt dir,side;
+  for (dir=0; dir<2; dir++) {
+    for (side=0; side<2; side++) {
+      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,0,0.0);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;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  if (dim < 1) {ierr = IGASetDim(iga,dim=2);CHKERRQ(ierr);}
+  ierr = IGASetUp(iga);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);
+
+  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 = SNESDestroy(&snes);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+
+  PetscBool flag = PETSC_FALSE;
+  PetscReal secs = -1;
+  ierr = PetscOptionsHasName(0,"-pause",&flag);CHKERRQ(ierr);
+  ierr = PetscOptionsGetReal(0,"-pause",&secs,0);CHKERRQ(ierr);
+  if (flag) {ierr = PetscSleep(secs);CHKERRQ(ierr);}
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File demo/BratuFJ.F90

+#include "petscconf.h"
+#if defined(PETSC_USE_COMPLEX)
+#define scalar complex
+#else
+#define scalar real
+#endif
+
+!module Bratu 
+!contains
+
+integer(kind=IGA_ERRCODE) &
+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)
+  do a = 1, p%nen
+     FF(a) = dot_product(grad_N(:,a),real(grad_u)) - N(a) * 1 * exp(real(u))
+  end do
+  ierr = 0
+end function Function
+
+integer(kind=IGA_ERRCODE) &
+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)
+  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))
+     end do
+  end do
+  ierr = 0
+end function Jacobian
+
+!end module Bratu 

File demo/makefile

           AdvectionDiffusion \
           Laplace \
           Poisson \
+          Bratu \
           PatternFormation \
           CahnHilliard2D \
           CahnHilliard3D \
 Poisson: Poisson.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<
+Bratu: Bratu.o BratuFJ.o chkopts
+	${CLINKER} -o $@ $< BratuFJ.o ${PETIGA_LIB}
+	${RM} -f $< BratuFJ.o
 AdvectionDiffusion: AdvectionDiffusion.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}
 	${RM} -f $<
 	-@${MPIEXEC} -n 1 ./PatternFormation ${OPTS} -ts_max_steps 2 -implicit
 runex5b_4:
 	-@${MPIEXEC} -n 4 ./PatternFormation ${OPTS} -ts_max_steps 2 -implicit
+runex6a_1:
+	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 1
+runex6a_2:
+	-@${MPIEXEC} -n 2 ./Bratu ${OPTS} -iga_dim 1
+runex6b_1:
+	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2
+runex6b_4:
+	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 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
 CahnHilliard2D := CahnHilliard2D.PETSc runex4_1 runex4_4 CahnHilliard2D.rm
 PatternFormation := PatternFormation.PETSc runex5a_1 runex5a_4 runex5b_1 runex5b_4 PatternFormation.rm
 
 Poisson      := $(Poisson1D) $(Poisson2D) $(Poisson3D)
 CahnHilliard := $(CahnHilliard2D) $(CahnHilliard3D)
 
-TESTEXAMPLES_C := $(InputOutput) $(L2Proj) $(Poisson) $(CahnHilliard) $(PatternFormation)
-TESTEXAMPLES_F :=
+TESTEXAMPLES_C := $(InputOutput) $(L2Proj) $(Poisson) $(Bratu) $(CahnHilliard) $(PatternFormation)
+TESTEXAMPLES_F := 
 TESTEXAMPLES_FORTRAN:=$(TESTEXAMPLES_F)
 testexamples:
 	-@${OMAKE} tree ACTION=testexamples_C PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR}

File include/petiga.h

 struct _n_IGAElement {
   PetscInt refct;
   /**/
-  PetscInt count;
-  PetscInt index;
   PetscInt start[3];
   PetscInt width[3];
   PetscInt ID[3];
-
+  /**/
+  PetscInt count;
+  PetscInt index;
+  /**/
   PetscInt nqp;
   PetscInt nen;
   PetscInt dof;
   /**/
   PetscInt count;
   PetscInt index;
-
+  /**/
   PetscInt nen;
   PetscInt dof;
   PetscInt dim;

File src/makefile

 #FFLAGS = ${FFLAGS_STD} -pedantic -Wall -Wextra -fcheck=all
 
 SOURCEH  = ../include/petiga.h petigapc.h petigagrid.h petigapart.h
-SOURCEC  = petiga.c petigareg.c petigaaxis.c petigarule.c petigabasis.c petigabound.c petigaelem.c petigapoint.c petigavec.c petigamat.c petigascl.c petigapcb.c petigapce.c petigaksp.c petigasnes.c petigats.c petigaio.c petigagrid.c petigapart.c
-SOURCEF1 = petigaftn.F90
-SOURCEF2 = petigabsp.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90 petigaint.f90
+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
+SOURCEF1 = petigaftn.F90 petigaint.F90
+SOURCEF2 = petigabsp.f90 petiga1d.f90 petiga2d.f90 petiga3d.f90
 SOURCEF  = ${SOURCEF1} ${SOURCEF2}
 OBJSC    = ${SOURCEC:.c=.o}
 OBJSF    = ${SOURCEF1:.F90=.o} ${SOURCEF2:.f90=.o}

File src/petigaftn.F90

 
 #include "petscconf.h"
 
-! PetscInt
+#define C_PETSC_ERRORCODE C_INT
+
 #if defined(PETSC_USE_64BIT_INDICES)
 #define C_PETSC_INT C_LONG_LONG
 #else
 #define C_PETSC_INT C_INT
 #endif
 
-! PetscReal
 #if defined(PETSC_USE_REAL_SINGLE)
 #define C_PETSC_REAL    C_FLOAT
 #define C_PETSC_COMPLEX C_FLOAT_COMPLEX
 #define C_PETSC_COMPLEX C_FLOAT128_COMPLEX
 #endif
 
-! PetscScalar
 #if defined(PETSC_USE_COMPLEX)
 #define C_PETSC_SCALAR  C_PETSC_COMPLEX
 #else
 #define C_PETSC_SCALAR  C_PETSC_REAL
 #endif
 
+#if defined(PETSC_USE_COMPLEX)
+#define scalar complex
+#else
+#define scalar real
+#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
   implicit none
+
+  !type, bind(C) :: IGA
+  !   private
+  !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
+  !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
+  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
+
+     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
+
+     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
+
+     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
+
+     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 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 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 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 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 IGAPointFormDer2
+     module procedure IGAPointFormHess_V
+     module procedure IGAPointFormHess_S
+  end interface IGAPointFormDer2
+
+  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
+
+  contains
+
+    integer(kind=IGA_ERRCODE) &
+    function IGAPointFormValue_S(p,U,v) bind(C) result (ierr) 
+      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
+
+    integer(kind=IGA_ERRCODE) &
+    function IGAPointFormValue_V(p,U,v) bind(C) result (ierr)
+      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
+
+    integer(kind=IGA_ERRCODE) &
+    function IGAPointFormGrad_S(p,U,v) bind(C) result (ierr)
+      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
+
+    integer(kind=IGA_ERRCODE) &
+    function IGAPointFormGrad_V(p,U,v) bind(C) result (ierr)
+      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
+
+    integer(kind=IGA_ERRCODE) &
+    function IGAPointFormHess_S(p,U,v) bind(C) result (ierr)
+      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
+
+    integer(kind=IGA_ERRCODE) &
+    function IGAPointFormHess_V(p,U,v) bind(C) result (ierr)
+      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
+
+    integer(kind=IGA_ERRCODE) &
+    function IGAPointFormDer3_S(p,U,v) bind(C) result (ierr)
+      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
+
+    integer(kind=IGA_ERRCODE) &
+    function IGAPointFormDer3_V(p,U,v) bind(C) result (ierr)
+      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
+
 end module PetIGA

File 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

File src/petigaint.F90

+! -*- f90 -*-
+
+#include "petscconf.h"
+#if defined(PETSC_USE_COMPLEX)
+#define scalar complex
+#else
+#define scalar real
+#endif
+
+subroutine IGA_GetPoint(nen,dim,N,C,X) &
+  bind(C, name="IGA_GetPoint")
+  use PetIGA
+  implicit none
+  integer(kind=IGA_INT ), intent(in),value :: nen,dim
+  real   (kind=IGA_REAL), intent(in)       :: N(nen)
+  real   (kind=IGA_REAL), intent(in)       :: C(dim,nen)
+  real   (kind=IGA_REAL), intent(out)      :: X(dim)
+  integer(kind=IGA_INT )  :: a
+  ! X = matmul(C,N)
+  X = 0
+  do a = 1, nen
+     X = X + N(a) * C(:,a)
+  end do
+end subroutine IGA_GetPoint
+
+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
+  ! V = matmul(N,transpose(U))
+  V = 0
+  do a = 1, nen
+     V = V + N(a) * U(:,a)
+  end do
+end subroutine IGA_GetValue
+
+subroutine IGA_GetGrad(nen,dof,dim,N,U,V) &
+  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
+  ! V = matmul(N,transpose(U))
+  V = 0
+  do a = 1, nen
+     do c = 1, dof
+        V(:,c) = V(:,c) + N(:,a) * U(c,a)
+     end do
+  end do
+end subroutine IGA_GetGrad
+
+subroutine IGA_GetHess(nen,dof,dim,N,U,V) &
+  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
+  ! V = matmul(N,transpose(U))
+  V = 0
+  do a = 1, nen
+     do i = 1, dof
+        V(:,i) = V(:,i) + N(:,a) * U(i,a)
+     end do
+  end do
+end subroutine IGA_GetHess
+
+subroutine IGA_GetDel2(nen,dof,dim,N,U,V) &
+  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
+  V = 0
+  do a = 1, nen
+     do c = 1, dof
+        do i = 1, dim
+           V(c) = V(c) + N(i,i,a) * U(c,a)
+        end do
+     end do
+  end do
+end subroutine IGA_GetDel2
+
+subroutine IGA_GetDer3(nen,dof,dim,N,U,V) &
+     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
+  ! V = matmul(N,transpose(U))
+  V = 0
+  do a = 1, nen
+     do i = 1, dof
+        V(:,i) = V(:,i) + N(:,a) * U(i,a)
+     end do
+  end do
+end subroutine IGA_GetDer3
+
+!subroutine IGA_GetDerivative(nen,dof,dim,der,N,U,V) &
+!  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
+!  ! V = matmul(N,transpose(U))
+!  V = 0
+!  do a = 1, nen
+!     do i = 1, dof
+!        V(:,i) = V(:,i) + N(:,a) * U(i,a)
+!     end do
+!  end do
+!end subroutine IGA_GetDerivative
+
+subroutine IGA_Interpolate(nen,dof,dim,der,N,U,V) &
+  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
+  ! V = matmul(N,transpose(U))
+  V = 0
+  do a = 1, nen
+     do i = 1, dof
+        V(:,i) = V(:,i) + N(:,a) * U(i,a)
+     end do
+  end do
+end subroutine IGA_Interpolate

File src/petigaint.f90

-! -*- f90 -*-
-
-subroutine IGA_GetPoint(nen,dim,N,C,X) &
-  bind(C, name="IGA_GetPoint")
-  use PetIGA
-  implicit none
-  integer(kind=IGA_INT ), intent(in),value :: nen,dim
-  real   (kind=IGA_REAL), intent(in)       :: N(nen)
-  real   (kind=IGA_REAL), intent(in)       :: C(dim,nen)
-  real   (kind=IGA_REAL), intent(out)      :: X(dim)
-  integer(kind=IGA_INT )  :: a
-  ! X = matmul(C,N)
-  X = 0
-  do a = 1, nen
-     X = X + N(a) * C(:,a)
-  end do
-end subroutine IGA_GetPoint
-
-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)
-  real   (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  real   (kind=IGA_SCALAR), intent(out)      :: V(dof)
-  integer(kind=IGA_INT   )  :: a, i
-  ! V = matmul(N,transpose(U))
-  V = 0
-  do a = 1, nen
-     V = V + N(a) * U(:,a)
-  end do
-end subroutine IGA_GetValue
-
-subroutine IGA_GetGrad(nen,dof,dim,N,U,V) &
-  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)
-  real   (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  real   (kind=IGA_SCALAR), intent(out)      :: V(dim,dof)
-  integer(kind=IGA_INT   )  :: a, c
-  ! V = matmul(N,transpose(U))
-  V = 0
-  do a = 1, nen
-     do c = 1, dof
-        V(:,c) = V(:,c) + N(:,a) * U(c,a)
-     end do
-  end do
-end subroutine IGA_GetGrad
-
-subroutine IGA_GetHess(nen,dof,dim,N,U,V) &
-  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)
-  real   (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  real   (kind=IGA_SCALAR), intent(out)      :: V(dim*dim,dof)
-  integer(kind=IGA_INT   )  :: a, i
-  ! V = matmul(N,transpose(U))
-  V = 0
-  do a = 1, nen
-     do i = 1, dof
-        V(:,i) = V(:,i) + N(:,a) * U(i,a)
-     end do
-  end do
-end subroutine IGA_GetHess
-
-subroutine IGA_GetDel2(nen,dof,dim,N,U,V) &
-  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)
-  real   (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  real   (kind=IGA_SCALAR), intent(out)      :: V(dof)
-  integer(kind=IGA_INT   )  :: a, c, i
-  V = 0
-  do a = 1, nen
-     do c = 1, dof
-        do i = 1, dim
-           V(c) = V(c) + N(i,i,a) * U(c,a)
-        end do
-     end do
-  end do
-end subroutine IGA_GetDel2
-
-subroutine IGA_GetDer3(nen,dof,dim,N,U,V) &
-     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)
-  real   (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  real   (kind=IGA_SCALAR), intent(out)      :: V(dim*dim*dim,dof)
-  integer(kind=IGA_INT   )  :: a, i
-  ! V = matmul(N,transpose(U))
-  V = 0
-  do a = 1, nen
-     do i = 1, dof
-        V(:,i) = V(:,i) + N(:,a) * U(i,a)
-     end do
-  end do
-end subroutine IGA_GetDer3
-
-!subroutine IGA_GetDerivative(nen,dof,dim,der,N,U,V) &
-!  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)
-!  real   (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-!  real   (kind=IGA_SCALAR), intent(out)      :: V(dim**der,dof)
-!  integer(kind=IGA_INT   )  :: a, i
-!  ! V = matmul(N,transpose(U))
-!  V = 0
-!  do a = 1, nen
-!     do i = 1, dof
-!        V(:,i) = V(:,i) + N(:,a) * U(i,a)
-!     end do
-!  end do
-!end subroutine IGA_GetDerivative
-
-subroutine IGA_Interpolate(nen,dof,dim,der,N,U,V) &
-  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)
-  real   (kind=IGA_SCALAR), intent(in)       :: U(dof,nen)
-  real   (kind=IGA_SCALAR), intent(out)      :: V(dim**der,dof)
-  integer(kind=IGA_INT   )  :: a, i
-  ! V = matmul(N,transpose(U))
-  V = 0
-  do a = 1, nen
-     do i = 1, dof
-        V(:,i) = V(:,i) + N(:,a) * U(i,a)
-     end do
-  end do
-end subroutine IGA_Interpolate