Commits

Huda Ibeid committed 5721087 Merge

Merged dalcinl/petiga into default

Comments (0)

Files changed (10)

   $ make test
 
 
+Scripting Support
+--------------
+
+PetIGA is designed to be efficient and as such, we do not directly do
+things like output VTK files suitable for viewing the solution. We do
+have routines which output the discretization information and solution
+vectors, but these are in a binary format to minimize I/O time. We
+have written a python package, `igakit
+<https://bitbucket.org/dalcinl/igakit>`_ which handles post-processing
+for visualization as well as geometry generation.
+
+
+Citation
+------
+
+If you find PetIGA helpful in conducting research projects, we would
+appreciate a citation to the following article::
+
+  @Article{PetIGA,
+    author = 	 {N. Collier, L. Dalcin, V.M. Calo},
+    title = 	 {{PetIGA}: High-Performance Isogeometric Analysis},
+    journal = 	 {arxiv},
+    year = 	 {2013},
+    number = 	 {1305.4452},
+    note = 	 {http://arxiv.org/abs/1305.4452}
+  }
+
+
 Acknowledgments
 ---------------
 
 PETIGA_LIB_DIR = ${PETIGA_DIR}/${PETSC_ARCH}/lib
 PETIGA_LIB     = ${CC_LINKER_SLFLAG}${PETIGA_LIB_DIR} -L${PETIGA_LIB_DIR} -lpetiga ${PETSC_LIB}
 
-CC_INCLUDES = ${PETSC_INCLUDE} ${PETIGA_INCLUDE}
-FC_INCLUDES = ${PETSC_INCLUDE} ${PETIGA_INCLUDE}
-CCPPFLAGS   = ${PETSC_CCPPFLAGS} ${PETIGA_INCLUDE}
-FCPPFLAGS   = ${PETSC_FCPPFLAGS} ${PETIGA_INCLUDE}
+CCPPFLAGS = ${PETSC_CCPPFLAGS} ${PETIGA_INCLUDE}
+FCPPFLAGS = ${PETSC_FCPPFLAGS} ${PETIGA_INCLUDE}
 
 INSTALL_LIB_DIR	= ${PETIGA_LIB_DIR}
-LIBNAME  = ${INSTALL_LIB_DIR}/${LIBBASE}.${AR_LIB_SUFFIX}
-SHLIBS   = libpetiga
+LIBNAME = ${INSTALL_LIB_DIR}/${LIBBASE}.${AR_LIB_SUFFIX}
+SHLIBS  = libpetiga

demo/TwoPhaseTwoComponent.c

   return;
 }
 
-#undef __FUNCT__
-#define __FUNCT__ "Permeability"
-PetscScalar Permeability(PetscScalar *p)
-{
-  if(p[0] < 100 && p[1] < 10){
-    return 1e-16;
-  }else if(p[0] >= 100 && p[1] < 10){
-    return 2e-16;
-  }else if(p[0] < 100 && p[1] >= 10){
-    return 3e-16;
-  }
-  return 4e-16;
-}
-
 #undef  __FUNCT__
 #define __FUNCT__ "Residual"
 PetscErrorCode Residual(IGAPoint p,PetscReal dt,
   PetscScalar Mh  = user->Mh;
   PetscScalar Mw  = user->Mw;
   PetscScalar Dlh = user->Dlh;
-  PetscScalar k   = user->k; //Permeability(p->point);
+  PetscScalar k   = user->k; 
 
   PetscInt a,i,nen,dim;
   IGAPointGetSizes(p,&nen,0,0);
 #define __FUNCT__ "GeometricAdaptivity"
 PetscErrorCode GeometricAdaptivity(TS ts,PetscReal t,Vec X,Vec Xdot, PetscReal *nextdt,PetscBool *ok,void *ctx)
 {
-  PetscErrorCode ierr;
-  PetscReal      dt;
-  ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
-  *ok = PETSC_TRUE;
-  *nextdt = PetscMin(1.2*dt,833.0);
-  return 0;
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "GeometricAdaptivity3"
-PetscErrorCode GeometricAdaptivity3(TS ts,PetscReal t,Vec X,Vec Xdot, PetscReal *nextdt,PetscBool *ok,void *ctx)
-{
   PetscErrorCode       ierr;
   PetscReal            dt;
   SNES                 snes;
   return 0;
 }
 
-
 #undef __FUNCT__
-#define __FUNCT__ "ICTest1"
-PetscErrorCode ICTest1(IGA iga,Vec U,AppCtx *user)
-{
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-
-  PetscScalar Mh = user->Mh;
-  PetscScalar H  = user->H;
-
-  DM da;
-  ierr = IGACreateNodeDM(iga,2,&da);CHKERRQ(ierr);
-  Field **u;
-  ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
-  DMDALocalInfo info;
-  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
-
-  PetscInt i,j;
-  for(i=info.xs;i<info.xs+info.xm;i++){
-    PetscReal x = (PetscReal)i / ( (PetscReal)(info.mx-1) );
-    for(j=info.ys;j<info.ys+info.ym;j++){
-      u[j][i].Pl = 10.;
-      if(x <= 0.5){
-	u[j][i].rholh = 15*Mh*H; // Pg = rholh/Mh
-      }else{
-	u[j][i].rholh = 25*Mh*H;
-      }
-    }
-  }
-  ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 
-  ierr = DMDestroy(&da);;CHKERRQ(ierr); 
-  PetscFunctionReturn(0); 
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "Test1"
-PetscErrorCode Test1()
-{
-
-  PetscErrorCode  ierr;
-  AppCtx user;
-  user.rholw    = 1000;
-  user.porosity = 0.3;
-  user.Pr       = 20;
-  user.n        = 1.54;
-  user.Slr      = 0.01;
-  user.Sgr      = 0.;
-  user.T        = 303;
-  user.H        = 0.765; 
-  user.Mw       = 1e-2;
-  user.Mh       = 2e-3;
-  user.mul      = 1e-8;
-  user.mug      = 9e-11;
-  user.Dlh      = 3e-9;
-  user.k        = 1e-16; 
-  PetscInt dim = 2, p = 1, N = 500,L = 1;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"r_","2c2p Options","2c2p");CHKERRQ(ierr);
-  ierr = PetscOptionsEnd();CHKERRQ(ierr);
-
-  IGA         iga;
-  IGAAxis     axis;
-  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
-  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
-  ierr = IGASetDof(iga,2);CHKERRQ(ierr);
-
-  ierr = IGAGetAxis(iga,0,&axis);CHKERRQ(ierr);
-  ierr = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis,N,0.0,L,0);CHKERRQ(ierr);
-  ierr = IGAGetAxis(iga,1,&axis);CHKERRQ(ierr);
-  ierr = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis,N/10,0.0,0.1*L,0);CHKERRQ(ierr);
-  
-  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-  ierr = IGASetUp(iga);CHKERRQ(ierr);
-  
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
-
-  TS     ts;
-  ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
-  ierr = TSSetDuration(ts,1000000,1.0e6);CHKERRQ(ierr);
-  ierr = TSSetTimeStep(ts,0.17);CHKERRQ(ierr);
-  ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
-  ierr = TSAlphaSetRadius(ts,0.75);CHKERRQ(ierr);
-  ierr = TSAlphaSetAdapt(ts,GeometricAdaptivity,NULL);CHKERRQ(ierr); 
-  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
-
-  Vec       U;
-  ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
-  ierr = ICTest1(iga,U,&user);CHKERRQ(ierr);
-  ierr = TSSolve(ts,U);CHKERRQ(ierr);
-
-  ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr);
-  ierr = IGAWriteVec(iga,U,"sol.dat");CHKERRQ(ierr);
-
-  ierr = VecDestroy(&U);CHKERRQ(ierr);
-  ierr = TSDestroy(&ts);CHKERRQ(ierr);
-  ierr = IGADestroy(&iga);CHKERRQ(ierr);
-
-  return 0;
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "ICTest3"
-PetscErrorCode ICTest3(IGA iga,Vec U,AppCtx *user)
+#define __FUNCT__ "InitialCondition"
+PetscErrorCode InitialCondition(IGA iga,Vec U,AppCtx *user)
 {
   PetscErrorCode ierr;
   PetscFunctionBegin;
 }
 
 #undef __FUNCT__
-#define __FUNCT__ "Test3Monitor"
-PetscErrorCode Test3Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *mctx)
+#define __FUNCT__ "Monitor"
+PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *mctx)
 {
   PetscErrorCode ierr;
   PetscFunctionBegin;
   AppCtx *user = (AppCtx *)mctx; 
-  IGA iga      = user->iga;
-  
-  // Painful crap needed to interpolate at the left-middle of the domain
-  IGAElement ele;
-  IGAPoint   pnt; 
+  IGA     iga  = user->iga;
+
+  // Pl,Pg,Sg computed at left middle  
   PetscScalar point[3] = {0,10,0};
-  ierr = IGAPointCreate(&pnt);CHKERRQ(ierr);
-  pnt->point = &(point[0]);
-  ierr = IGAElementCreate(&ele);CHKERRQ(ierr);
-  ierr = IGAElementInit(ele,iga);CHKERRQ(ierr);
-  if(IGALocateElement(iga,pnt->point,ele)){
-    ierr = IGAPointInit(pnt,ele);CHKERRQ(ierr);
-    ierr = IGAPointEval(iga,pnt);CHKERRQ(ierr);
-  }else{
-    SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Point not on IGA.");
-  }
-  Vec                localU;
-  const PetscScalar *arrayU;
-  PetscScalar       *Ul;
-  ierr = IGAGetLocalVecArray(iga,U,&localU,&arrayU);CHKERRQ(ierr);
-  ierr = IGAElementGetWorkVal(ele,&Ul);CHKERRQ(ierr);
-  ierr = IGAElementGetValues(ele,arrayU,Ul);CHKERRQ(ierr);
-
-  // Pl,Pg,Sg computed at left middle
   PetscScalar sol[2];
-  ierr = IGAPointFormValue(pnt,Ul,&sol[0]);
+  ierr = IGAInterpolate(iga,U,point,&sol[0],NULL);CHKERRQ(ierr);
   PetscScalar Pl,Pg,Sg;
   Pl = sol[0];
   Pg = sol[1]/(user->Mh*user->H);
 
   PetscReal dt;
   TSGetTimeStep(ts,&dt);
+  if(step == 0){
+    PetscPrintf(PETSC_COMM_WORLD,"#%11s %12s %12s %12s %12s %12s %12s\n","Time","dt","Pl(xL)","Pg(xL)","Sg(xL)","fluxw(xR)","fluxh(xR)");
+  }
   PetscPrintf(PETSC_COMM_WORLD,"%.6e %.6e %.6e %.6e %.6e %.6e %.6e\n",t,dt,Pl,Pg,Sg,fluxw,fluxh);
-  
-  ierr = IGARestoreLocalVecArray(iga,U,&localU,&arrayU);CHKERRQ(ierr);
-  ierr = IGAElementDestroy(&ele);CHKERRQ(ierr);
-  ierr = IGAPointDestroy(&pnt);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
 
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
 
-#undef __FUNCT__
-#define __FUNCT__ "Test3"
-PetscErrorCode Test3()
-{
   PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
   AppCtx user;
   user.rholw    = 1000;
   ierr = TSSetDuration(ts,1000000,1.0e6);CHKERRQ(ierr);
   ierr = TSSetTimeStep(ts,10.0);CHKERRQ(ierr);
   ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
-  ierr = TSAlphaSetRadius(ts,0.75);CHKERRQ(ierr);
-  ierr = TSAlphaSetAdapt(ts,GeometricAdaptivity3,NULL);CHKERRQ(ierr); 
-  ierr = TSMonitorSet(ts,Test3Monitor,&user,NULL);CHKERRQ(ierr);
+  ierr = TSAlphaSetRadius(ts,0.95);CHKERRQ(ierr);
+  ierr = TSAlphaSetAdapt(ts,GeometricAdaptivity,NULL);CHKERRQ(ierr); 
+  ierr = TSMonitorSet(ts,Monitor,&user,NULL);CHKERRQ(ierr);
   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
 
   Vec       U;
   ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
-  ierr = ICTest3(iga,U,&user);CHKERRQ(ierr);
+  ierr = InitialCondition(iga,U,&user);CHKERRQ(ierr);
   ierr = TSSolve(ts,U);CHKERRQ(ierr);
 
-  ierr = IGAWrite(iga,"iga_test3.dat");CHKERRQ(ierr);
-  ierr = IGAWriteVec(iga,U,"ss_test3.dat");CHKERRQ(ierr);
+  ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr);
+  ierr = IGAWriteVec(iga,U,"ss.dat");CHKERRQ(ierr);
 
   ierr = VecDestroy(&U);CHKERRQ(ierr);
   ierr = TSDestroy(&ts);CHKERRQ(ierr);
   ierr = IGADestroy(&iga);CHKERRQ(ierr);
 
-  return ierr;
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "main"
-int main(int argc, char *argv[]) {
-
-  PetscErrorCode  ierr;
-  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
-
-  //Test1();
-  //Test2();
-  ierr = Test3();CHKERRQ(ierr);
-
   ierr = PetscFinalize();CHKERRQ(ierr);
   return 0;
 }
 PETSC_EXTERN PetscErrorCode IGASetUserIEFunction(IGA iga,IGAUserIEFunction IEFunction,void *ctx);
 PETSC_EXTERN PetscErrorCode IGASetUserIEJacobian(IGA iga,IGAUserIEJacobian IEJacobian,void *ctx);
 
+PETSC_EXTERN PetscBool IGALocateElement(IGA iga,PetscReal *pnt,IGAElement element);
+PETSC_EXTERN PetscErrorCode IGAPointEval(IGA iga,IGAPoint point);
+PETSC_EXTERN PetscErrorCode IGAInterpolate(IGA iga,Vec U,PetscReal p[],PetscScalar u[],PetscScalar du[]);
+
 #undef  DMIGA
 #define DMIGA "iga"
 PETSC_EXTERN PetscErrorCode IGACreateWrapperDM(IGA iga,DM *dm);
 
 /* ---------------------------------------------------------------- */
 
-PETSC_EXTERN PetscBool IGALocateElement(IGA iga,PetscReal *pnt,IGAElement element);
-PETSC_EXTERN PetscErrorCode IGAPointEval(IGA iga,IGAPoint point);
-
-/* ---------------------------------------------------------------- */
-
 #if defined(PETSC_USE_DEBUG)
 #define IGACheckSetUp(iga,arg) do {                                      \
     if (PetscUnlikely(!(iga)->setup))                                    \
 cmake_cc_path =-DCMAKE_C_COMPILER:FILEPATH=${CC}
 cmake_cc_flags=-DCMAKE_C_FLAGS:STRING='${PCC_FLAGS} ${CFLAGS} ${PETSCFLAGS} ${CPP_FLAGS} ${CPPFLAGS}'
 endif
+ifneq (${FC},)
 cmake_fc_path =-DCMAKE_Fortran_COMPILER:FILEPATH=${FC}
 cmake_fc_flags=-DCMAKE_Fortran_FLAGS:STRING='${FC_FLAGS} ${FFLAGS} ${PETSCFLAGS} ${FPP_FLAGS} ${FPPFLAGS}'
+endif
 cmake_cc=${cmake_cc_path} ${cmake_cc_flags} ${cmake_cc_clang}
 cmake_fc=${cmake_fc_path} ${cmake_fc_flags}
-${PETIGA_DIR}/${PETSC_ARCH}/CMakeCache.txt: CMakeLists.txt ${PETIGA_DIR}/${PETSC_ARCH}/conf
+${PETIGA_DIR}/${PETSC_ARCH}/CMakeCache.txt: CMakeLists.txt
 	-@${RM} -r ${PETIGA_DIR}/${PETSC_ARCH}/CMakeCache.txt
 	-@${RM} -r ${PETIGA_DIR}/${PETSC_ARCH}/CMakeFiles
 	-@${RM} -r ${PETIGA_DIR}/${PETSC_ARCH}/Makefile
 	-@${RM} -r ${PETIGA_DIR}/${PETSC_ARCH}/cmake_install.cmake
-	@cd ${PETIGA_DIR}/${PETSC_ARCH} && ${CMAKE} ${PETIGA_DIR} ${cmake_cc} ${cmake_fc} 2>&1 > ${PETIGA_DIR}/${PETSC_ARCH}/conf/cmake.log
+	@${MKDIR} ${PETIGA_DIR}/${PETSC_ARCH}
+	@cd ${PETIGA_DIR}/${PETSC_ARCH} && ${CMAKE} ${PETIGA_DIR} ${cmake_cc} ${cmake_fc}
 cmake-boot: ${PETIGA_DIR}/${PETSC_ARCH}/CMakeCache.txt
 cmake-build: cmake-boot
-	@cd ${PETIGA_DIR}/${PETSC_ARCH} && ${OMAKE} -j ${MAKE_NP} 2>&1
+	@cd ${PETIGA_DIR}/${PETSC_ARCH} && ${OMAKE} -j ${MAKE_NP}
 cmake-clean:
 	-@${RM} -r ${PETIGA_DIR}/${PETSC_ARCH}/CMakeCache.txt
 	-@${RM} -r ${PETIGA_DIR}/${PETSC_ARCH}/CMakeFiles
 	-@echo "============================================="
 	@${OMAKE} cmake-build PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR} 2>&1 | tee ./${PETSC_ARCH}/conf/make.log
 	-@echo "============================================="
-.PHONY: cmake-boot cmake-build all-cmake
+.PHONY: cmake-boot cmake-build cmake-clean all-cmake
 
 #
 # Check if PETSC_DIR variable specified is valid
   #elif defined(PETSC_USE_REAL_DOUBLE)
     #define LAPACKgetri_ zgetri_
   #else /* (PETSC_USE_REAL_QUAD) */
-    #error "LAPACKgetri_ not defined for quad complex"
+    #define LAPACKgetri_ wgetri_
   #endif
 #endif
 EXTERN_C_BEGIN
   end interface IGA_InvGradGeomMap
 
 
+  interface IGA_Basis0
+     module procedure IGA_Basis0
+  end interface IGA_Basis0
+
+  interface IGA_Basis1
+     module procedure IGA_Basis1
+  end interface IGA_Basis1
+
+  interface IGA_Basis2
+     module procedure IGA_Basis2
+  end interface IGA_Basis2
+
+  interface IGA_Basis3
+     module procedure IGA_Basis3
+  end interface IGA_Basis3
+
+
   interface IGA_Shape0
      module procedure IGA_Shape0
   end interface IGA_Shape0
       call c2f(p%normal,N,(/p%dim/))
     end function IGA_Normal
 
+    function IGA_Basis0(p) result(N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      real(kind=IGA_REAL_KIND), pointer :: N(:)
+      call c2f(p%basis(0),N,(/p%nen/))
+    end function IGA_Basis0
+
+    function IGA_Basis1(p) result(N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      real(kind=IGA_REAL_KIND), pointer :: N(:,:)
+      call c2f(p%basis(1),N,(/p%dim,p%nen/))
+    end function IGA_Basis1
+
+    function IGA_Basis2(p) result(N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      real(kind=IGA_REAL_KIND), pointer :: N(:,:,:)
+      call c2f(p%basis(2),N,(/p%dim,p%dim,p%nen/))
+    end function IGA_Basis2
+
+    function IGA_Basis3(p) result(N)
+      use ISO_C_BINDING, only: c2f => C_F_POINTER
+      implicit none
+      type(IGAPoint), intent(in) :: p
+      real(kind=IGA_REAL_KIND), pointer :: N(:,:,:,:)
+      call c2f(p%basis(3),N,(/p%dim,p%dim,p%dim,p%nen/))
+    end function IGA_Basis3
+
     function IGA_Shape0(p) result(N)
       use ISO_C_BINDING, only: c2f => C_F_POINTER
       implicit none
 
 .keywords: IGA, locate, point, element
 @*/
-PetscBool IGALocateElement(IGA iga,PetscScalar *pnt,IGAElement element)
+PetscBool IGALocateElement(IGA iga,PetscReal *pnt,IGAElement element)
 {
   PetscErrorCode ierr;
   PetscInt i,j,e,m,dim=iga->dim,*ID = element->ID;
 +  iga - the IGA context
 -  point - the IGAPoint context
 
-   Notes:
-   The point assumes you have already called IGALocateElement and it
-   returned true. This function is intended to be used when
-   integrating on a manifold embedded in the IGA domain.
+   Notes: The point assumes that the point is contained in an element
+   on this processors partition (call IGALocateElement)
 
    Level: devel
 
   for(i=0;i<element->dim;i++){ ierr = PetscFree(basis[i]);CHKERRQ(ierr); }
   PetscFunctionReturn(0);
 }
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAInterpolate"
+/*@
+   IGAInterpolate - Interpolates a vector given the IGA at the
+   requested point in the domain.
+
+   Input Parameters:
++  iga - the IGA context
+.  U   - the vector
+.  p   - the parametric location in the IGA [iga->dim]
+.  u   - pointer to the interpolated values [iga->dof]
+-  du  - pointer to the gradient of the interpolated values [iga->dof][iga->dim]
+
+   Notes: This function is intended to be used in a monitor function
+   if a few values of the solution are needed at specific
+   locations. Using it for interpolating a large number of points will
+   be slow as it needs to create and destroy several objects
+   internally.
+
+   Level: normal
+
+.keywords: interpolate
+@*/
+PetscErrorCode IGAInterpolate(IGA iga,Vec U,PetscReal p[],PetscScalar u[],PetscScalar du[])
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(U,VEC_CLASSID,2);
+  PetscValidRealPointer(p,3);
+  if(u)  PetscValidScalarPointer(u ,4);
+  if(du) PetscValidScalarPointer(du,5);
+  PetscErrorCode     ierr;
+  IGAElement         ele;
+  IGAPoint           pnt;
+  MPI_Comm           comm;
+  PetscMPIInt        rank,lroot=0,root=0;
+  PetscMPIInt        lfound=0,found=0;
+  Vec                localU;
+  const PetscScalar *arrayU;
+  PetscScalar       *Ul;
+
+  /* get comm and rank */
+  ierr = IGAGetComm(iga,&comm);CHKERRQ(ierr);
+  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
+
+  /* create the IGAPoint and IGAElement */
+  ierr = IGAPointCreate(&pnt);CHKERRQ(ierr);
+  ierr = IGAElementCreate(&ele);CHKERRQ(ierr);
+  ierr = IGAElementInit(ele,iga);CHKERRQ(ierr);
+  pnt->point = p;
+  pnt->dof   = iga->dof;
+  pnt->dim   = iga->dim;
+
+  /* if point is on this processor, this rank is root */
+  ierr = IGAGetLocalVecArray(iga,U,&localU,&arrayU);CHKERRQ(ierr);
+  if(IGALocateElement(iga,pnt->point,ele)){
+    lroot  = rank;
+    lfound = 1;
+    ierr = IGAPointInit(pnt,ele);CHKERRQ(ierr);
+    ierr = IGAPointEval(iga,pnt);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVal(ele,&Ul);CHKERRQ(ierr);
+    ierr = IGAElementGetValues(ele,arrayU,Ul);CHKERRQ(ierr);
+    if(u)  IGA_GetValue(pnt->nen,pnt->dof,pnt->shape[0],Ul,u);
+    if(du) IGA_GetGrad(pnt->nen,pnt->dof,pnt->dim,pnt->shape[1],Ul,du);
+  }
+  ierr = IGARestoreLocalVecArray(iga,U,&localU,&arrayU);CHKERRQ(ierr);
+  ierr = MPI_Allreduce(&lfound,&found,1,MPI_INT,MPI_SUM,comm);
+  if(found == 0){
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Point not found on any partition of the domain");
+  }
+  ierr = MPI_Allreduce(&lroot,&root,1,MPI_INT,MPI_SUM,comm);
+
+  /* broadcast result so all processors have same data */
+  if(u)  { ierr = MPI_Bcast( u,pnt->dof,         MPIU_SCALAR,root,comm);CHKERRQ(ierr); }
+  if(du) { ierr = MPI_Bcast(du,pnt->dof*pnt->dim,MPIU_SCALAR,root,comm);CHKERRQ(ierr); }
+  ierr = IGAElementDestroy(&ele);CHKERRQ(ierr);
+  ierr = IGAPointDestroy(&pnt);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+#include "petiga.h"
+#include "time.h"
+
+#define SQ(A) ((A)*(A))
+
+PetscScalar CrossTerm(PetscReal x, PetscReal y)
+{
+  return x*y;
+}
+
+typedef struct {
+  PetscScalar (*Function)(PetscReal x, PetscReal y);
+} AppCtx;
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal x = p->point[0];
+  PetscReal y = p->point[1];
+  PetscScalar f = user->Function(x,y);
+
+  PetscReal *N = p->shape[0];
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na = N[a];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb = N[b];
+      K[a*nen+b] = Na * Nb;
+    }
+    F[a] = Na * f;
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  AppCtx user;
+  user.Function = CrossTerm;
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,2);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGASetUserSystem(iga,System,&user);CHKERRQ(ierr);
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+
+  KSP ksp;
+  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+
+  PetscReal  pnt[2];
+  PetscScalar sol=0,gsol[2],err,gerr;
+  srand(time(NULL));
+  PetscInt i;
+  for(i=0;i<10;i++){
+    pnt[0] = rand();
+    pnt[1] = rand();
+    pnt[0] /= RAND_MAX;
+    pnt[1] /= RAND_MAX;
+    PetscPrintf(PETSC_COMM_WORLD,"Evaluating solution at x = (%f,%f)\n",pnt[0],pnt[1]);
+    ierr = IGAInterpolate(iga,x,pnt,&sol,gsol);
+    err = fabs(sol-user.Function(pnt[0],pnt[1]));
+    gerr = sqrt(SQ(gsol[0]-pnt[1])+SQ(gsol[1]-pnt[0]));
+    PetscPrintf(PETSC_COMM_WORLD,"\t     Solution  =  %e                Error = %e\n",sol,err);
+    PetscPrintf(PETSC_COMM_WORLD,"\tgrad(Solution) = [%e %e]  Error = %e\n",gsol[0],gsol[1],gerr);
+  }
+
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = MatDestroy(&A);CHKERRQ(ierr);
+  ierr = VecDestroy(&x);CHKERRQ(ierr);
+  ierr = VecDestroy(&b);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}
 
 OPTS=-nox -malloc_debug -malloc_dump
 
+Test_Eval: Test_Eval.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
 
 IGACreate: IGACreate.o chkopts
 	${CLINKER} -o $@ $< ${PETIGA_LIB}