Commits

Lisandro Dalcin committed 946f2e9

Huge changes, remove IGABoundary and add IGAForm, many renames

  • Participants
  • Parent commits 81bfb26

Comments (0)

Files changed (48)

File demo/AdaptiveL2Projection.c

   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,System,NULL);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
   PetscScalar scalar[2];
-  ierr  = IGAFormScalar(iga,x,2,&scalar[0],Error,NULL);CHKERRQ(ierr);
+  ierr  = IGAComputeScalar(iga,x,2,&scalar[0],Error,NULL);CHKERRQ(ierr);
   PetscReal errors[2];
   errors[0] = PetscSqrtReal(PetscRealPart(scalar[0]));
   errors[1] = PetscSqrtReal(PetscRealPart(scalar[1]));
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,System,NULL);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;

File demo/AdvectionDiffusion.c

   AppCtx user;
   ierr = ComputeWind(dim,Pe,dir,user.wind);CHKERRQ(ierr);
 
-  IGABoundary bnd;
   PetscReal ul=1.0, ur=0.0;
   for (i=0; i<dim; i++) {
     if (dir[i] == 0.0) continue;
-    ierr = IGAGetBoundary(iga,i,0,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetValue(bnd,0,ul);CHKERRQ(ierr);
-    ierr = IGAGetBoundary(iga,i,1,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetValue(bnd,0,ur);CHKERRQ(ierr);
+    ierr = IGASetBoundaryValue(iga,i,0,0,ul);CHKERRQ(ierr);
+    ierr = IGASetBoundaryValue(iga,i,1,0,ur);CHKERRQ(ierr);
   }
 
   Mat A;
   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 = IGASetFormSystem(iga,System,&user);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;

File demo/BoundaryIntegral.c

 #include "petiga.h"
 
 typedef struct {
-  PetscInt dir;
+  PetscInt axis;
   PetscInt side;
 } AppCtx;
 
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "System"
-PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+#define __FUNCT__ "Laplace"
+PetscErrorCode Laplace(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt dim = p->dim;
   PetscInt nen = p->nen;
   return 0;
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  if (!p->atboundary)
+    return Laplace(p,K,F,ctx);
+  else
+    return Neumann(p,K,F,ctx);
+}
+
+
 PETSC_STATIC_INLINE
 PetscReal DEL2(PetscInt dim,const PetscReal a[dim][dim])
 {
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "SystemCollocation"
-PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+#define __FUNCT__ "LaplaceCollocation"
+PetscErrorCode LaplaceCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   PetscInt nen = p->nen;
   PetscInt dim = p->dim;
 {
   PetscInt nen = p->nen;
   const PetscReal *N0 = (typeof(N0)) p->shape[0];
-
   PetscInt a;
   for (a=0; a<nen; a++)
     K[a] = N0[a];
 PetscErrorCode Neumann0Collocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {(void)NeumannCollocation(p,K,F,ctx); F[0] = 0.0; return 0;}
 
+
+#undef  __FUNCT__
+#define __FUNCT__ "SystemCollocation"
+PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  if (!p->atboundary) return LaplaceCollocation(p,K,F,ctx);
+
+  AppCtx *user = (AppCtx *)ctx;
+  PetscInt axis,side;
+  IGAPointAtBoundary(p,&axis,&side);
+  if (axis == user->axis) {
+    if (side == user->side)
+      return NeumannCollocation(p,K,F,ctx);
+    else
+      return DirichletCollocation(p,K,F,ctx);
+  } else
+    return Neumann0Collocation(p,K,F,ctx);
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "Error"
 PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
   IGAPointFormValue(p,U,&u);
   PetscReal x;
   if (user->side == 0)
-    x = 1 - p->point[user->dir] + 1;
+    x = 1 - p->point[user->axis] + 1;
   else
-    x = p->point[user->dir] + 1;
+    x = p->point[user->axis] + 1;
   PetscReal e = u - x;
   S[0] = e*e;
   return 0;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
   AppCtx user;
-  user.dir  = 0;
+  user.axis = 0;
   user.side = 1;
 
   PetscBool print_error = PETSC_FALSE;
   PetscBool check_error = PETSC_FALSE;
   PetscBool draw = PETSC_FALSE;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","BoundaryIntegral Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt ("-dir", "Neuman BC direction",__FILE__,user.dir, &user.dir, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt ("-axis","Neuman BC direction",__FILE__,user.axis,&user.axis,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsInt ("-side","Neuman BC side",     __FILE__,user.side,&user.side,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-print_error","Prints the L2 error of the solution",__FILE__,print_error,&print_error,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-check_error","Checks the L2 error of the solution",__FILE__,check_error,&check_error,NULL);CHKERRQ(ierr);
 
   PetscInt dim;
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
-  PetscInt d =  !user.side;
-  PetscInt n = !!user.side;
+
+  IGAForm form;
+  ierr = IGAGetForm(iga,&form);CHKERRQ(ierr);
   if (!iga->collocation) {
-    IGABoundary bnd;
-    ierr = IGAGetBoundary(iga,user.dir,d,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetValue(bnd,0,1.0);CHKERRQ(ierr);
-    ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetUserSystem(bnd,Neumann,NULL);CHKERRQ(ierr);
+    PetscInt d = !user.side;
+    PetscInt n = !d;
+    ierr = IGAFormSetSystem(form,System,&user);CHKERRQ(ierr);
+    ierr = IGAFormSetBoundaryValue(form,user.axis,d,0,1.0);CHKERRQ(ierr);
+    ierr = IGAFormSetBoundaryForm (form,user.axis,n,PETSC_TRUE);CHKERRQ(ierr);
   } else {
-    IGABoundary bnd;
-    PetscInt dir,side;
-    for (dir=0; dir<dim; dir++) {
-      for (side=0; side<2; side++) {
-        ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-        ierr = IGABoundarySetUserSystem(bnd,Neumann0Collocation,NULL);CHKERRQ(ierr);
-      }
-    }
-    ierr = IGAGetBoundary(iga,user.dir,d,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetUserSystem(bnd,DirichletCollocation,NULL);CHKERRQ(ierr);
-    ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetUserSystem(bnd,NeumannCollocation,NULL);CHKERRQ(ierr);
+    PetscInt i,s;
+    ierr = IGAFormSetSystem(form,SystemCollocation,&user);CHKERRQ(ierr);
+    for (i=0; i<dim; i++)
+      for (s=0; s<2; s++)
+        {ierr = IGAFormSetBoundaryForm(form,i,s,PETSC_TRUE);CHKERRQ(ierr);}
   }
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  if (!iga->collocation) {
-    ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
-  } else {
-    ierr = IGASetUserSystem(iga,SystemCollocation,NULL);CHKERRQ(ierr);
-  }
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
   PetscScalar error = 0;
-  ierr = IGAFormScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
 
   if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);}

File demo/Bratu.c

   PetscInt dim,dir,side;
   for (dir=0; dir<3; dir++) {
     for (side=0; side<2; side++) {
-      IGABoundary bnd;
       PetscInt    field = 0;
       PetscScalar value = 0.0;
-      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,field,value);CHKERRQ(ierr);
+      ierr = IGASetBoundaryValue(iga,dir,side,field,value);CHKERRQ(ierr);
     }
   }
   ierr = IGASetFromOptions(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);
+    ierr = IGASetFormFunction(iga,Bratu_Function,&ctx);CHKERRQ(ierr);
+    ierr = IGASetFormJacobian(iga,Bratu_Jacobian,&ctx);CHKERRQ(ierr);
   } else {
-    ierr = IGASetUserIFunction(iga,Bratu_IFunction,&ctx);CHKERRQ(ierr);
-    ierr = IGASetUserIJacobian(iga,Bratu_IJacobian,&ctx);CHKERRQ(ierr);
+    ierr = IGASetFormIFunction(iga,Bratu_IFunction,&ctx);CHKERRQ(ierr);
+    ierr = IGASetFormIJacobian(iga,Bratu_IJacobian,&ctx);CHKERRQ(ierr);
   }
 
   Vec x;

File demo/CahnHilliard2D.c

   AppCtx *user = (AppCtx *)mctx;
 
   PetscScalar stats[3] = {0,0,0};
-  ierr = IGAFormScalar(user->iga,U,3,&stats[0],Stats,mctx);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(user->iga,U,3,&stats[0],Stats,mctx);CHKERRQ(ierr);
 
   PetscReal dt;
   TSGetTimeStep(ts,&dt);
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
 
 
   TS ts;

File demo/CahnHilliard3D.c

   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
   th->time[0]=t-dt;
   PetscScalar energy;
-  ierr = IGAFormScalar(th->iga,X,1,&energy,Stats,ctx);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(th->iga,X,1,&energy,Stats,ctx);CHKERRQ(ierr);
   th->energy[0] = PetscRealPart(energy);
 
   if (!th->E){
 			 iga->elem_sizes[2]*iga->elem_sizes[2]);
   user.lambda = user.tau*h*h; /* mesh size parameter */
 
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
 
   TS ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);

File demo/ClassicalShell.c

   ierr = IGARead(iga,filename);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  IGABoundary bnd;
   // Boundary conditions on u = 0, v = [0:1]
-  ierr = IGAGetBoundary(iga,0,0,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,1,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,2,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,4,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,1,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,2,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,4,0.0);CHKERRQ(ierr);
   // Boundary conditions on u = 1, v = [0:1]
-  ierr = IGAGetBoundary(iga,0,1,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,3,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,1,0,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,1,3,0.0);CHKERRQ(ierr);
   // Boundary conditions on u = [0:1], v = 0
-  ierr = IGAGetBoundary(iga,1,0,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,1,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,4,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,1,0,1,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,1,0,4,0.0);CHKERRQ(ierr);
   // Boundary conditions on u = [0:1], v = 1
-  ierr = IGAGetBoundary(iga,1,1,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,1,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,4,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,1,1,1,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,1,1,4,0.0);CHKERRQ(ierr);
 
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   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 = IGASetFormSystem(iga,System,&user);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
   //MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
   //PetscBool symm=1;

File demo/ElasticRod.c

   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   /* Boundary conditions, U[left] = U[right] = 0 */
-  IGABoundary bnd;
   PetscScalar U_left = 0.0, U_right = 0.0;
   /* left,  (dir=0, side=0) */
-  ierr = IGAGetBoundary(iga,0,0,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,U_left);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,0,U_left);CHKERRQ(ierr);
   /* right  (dir=0, side=1) */
-  ierr = IGAGetBoundary(iga,0,1,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,U_right);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,1,0,U_right);CHKERRQ(ierr);
 
   /* Residual and Tangent user routines */
-  ierr = IGASetUserIFunction2(iga,ElasticRod_IFunction,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian2(iga,ElasticRod_IJacobian,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction2(iga,ElasticRod_IFunction,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian2(iga,ElasticRod_IJacobian,&user);CHKERRQ(ierr);
 
   /* Timestepper, t_final=5.0, delta_t = 0.01 */
   TS ts;

File demo/Elasticity3D.c

   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   // Set boundary conditions
-  IGABoundary bnd;
-  ierr = IGAGetBoundary(iga,0,0,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,1,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,2,0.0);CHKERRQ(ierr);
-  ierr = IGAGetBoundary(iga,0,1,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,1.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,0,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,1,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,2,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,1,0,1.0);CHKERRQ(ierr);
 
   // Create linear system
   Mat A;
   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 = IGASetFormSystem(iga,System,&user);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
   
   // Solve

File demo/HyperElasticity.c

   }
 
   // Set boundary conditions
-  IGABoundary bnd;
-  ierr = IGAGetBoundary(iga,0,0,&bnd);CHKERRQ(ierr); // u = [0,0,0] @ x = [0,:,:]
-  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,1,0.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,2,0.0);CHKERRQ(ierr);
-
-  ierr = IGAGetBoundary(iga,0,1,&bnd);CHKERRQ(ierr); // ux = 1 @ x = [1,:,:]
-  ierr = IGABoundarySetValue(bnd,0,0.1/((PetscScalar)nsteps));CHKERRQ(ierr);
+  //   u = [0,0,0] @ x = [0,:,:]
+  ierr = IGASetBoundaryValue(iga,0,0,0,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,1,0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,0,2,0.0);CHKERRQ(ierr);
+  //   ux = 1 @ x = [1,:,:]
+  ierr = IGASetBoundaryValue(iga,0,1,0,0.1/((PetscScalar)nsteps));CHKERRQ(ierr);
 
   // Setup the nonlinear solver
   SNES snes;
   Vec U,Utotal;
-  ierr = IGASetUserFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
+  ierr = IGASetFormFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
   ierr = IGACreateSNES(iga,&snes);CHKERRQ(ierr);
   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);

File demo/L2Proj1D.c

   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 = IGASetFormSystem(iga,System,&user);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
   PetscScalar error = 0;
-  ierr = IGAFormScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
   PetscBool print_error = PETSC_FALSE;
   ierr = PetscOptionsGetBool(0,"-error",&print_error,0);CHKERRQ(ierr);

File demo/L2Proj2D.c

   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 = IGASetFormSystem(iga,System,&user);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
   PetscScalar error = 0;
-  ierr = IGAFormScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
   PetscBool print_error = PETSC_FALSE;
   ierr = PetscOptionsGetBool(0,"-error",&print_error,0);CHKERRQ(ierr);

File demo/L2Projection.c

   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,System,NULL);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
   PetscScalar error = 0;
-  ierr = IGAFormScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
   PetscBool print_error = PETSC_FALSE;
   ierr = PetscOptionsGetBool(0,"-print_error",&print_error,0);CHKERRQ(ierr);

File demo/Laplace.c

   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
   if (collocation) {ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr);}
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-  if (iga->dim < 1) {ierr = IGASetDim(iga,2);CHKERRQ(ierr);}
+  if (iga->dim<1) {ierr = IGASetDim(iga,2);CHKERRQ(ierr);}
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   // Set boundary conditions
 
-  PetscInt dim,dir;
+  PetscInt dim,dir,side;
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   for (dir=0; dir<dim; dir++) {
-    IGABoundary bnd;
     PetscScalar value = 1.0;
     PetscScalar load  = 0.0;
-    ierr = IGAGetBoundary(iga,dir,0,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
-    ierr = IGAGetBoundary(iga,dir,1,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetLoad(bnd,0,load);CHKERRQ(ierr);
+    ierr = IGASetBoundaryValue(iga,dir,side=0,0,value);CHKERRQ(ierr);
+    ierr = IGASetBoundaryLoad (iga,dir,side=1,0,load );CHKERRQ(ierr);
   }
 
   // Assemble
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   if (!iga->collocation) {
-    ierr = IGASetUserSystem(iga,SystemGalerkin,NULL);CHKERRQ(ierr);
+    ierr = IGASetFormSystem(iga,SystemGalerkin,NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
   } else {
-    ierr = IGASetUserSystem(iga,SystemCollocation,NULL);CHKERRQ(ierr);
+    ierr = IGASetFormSystem(iga,SystemCollocation,NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_FALSE);CHKERRQ(ierr);
   }
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
   // Various post-processing options
 
   PetscScalar error = 0;
-  ierr = IGAFormScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
 
   if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);}

File demo/LoggChallenge.c

   ierr = IGASetDim(iga,2);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);
-    }
-  }
+  for (dir=0; dir<2; dir++)
+    for (side=0; side<2; side++)
+      {ierr = IGASetBoundaryValue(iga,dir,side,0,0.0);CHKERRQ(ierr);}
 
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,System,NULL);CHKERRQ(ierr);
 
   PetscLogDouble ta1,ta2,ta;
   ierr = PetscTime(&ta1);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   PetscScalar error = 0;
-  ierr = IGAFormScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
 
   ierr = PetscPrintf(PETSC_COMM_WORLD,

File demo/NavierStokesKorteweg2D.c

   AppCtx *user = (AppCtx *)mctx;
 
   PetscScalar scalar = 0.;
-  ierr = IGAFormScalar(user->iga,U,1,&scalar,Energy,mctx);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(user->iga,U,1,&scalar,Energy,mctx);CHKERRQ(ierr);
   PetscReal energy = PetscRealPart(scalar);
 
   PetscReal dt;
   ierr = IGAWrite(iga,"iga.dat");CHKERRQ(ierr);
   user.iga = iga;
 
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
 
   TS ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);

File demo/NavierStokesVMS.c

   //ierr = IGAAxisInitUniform(axis2,N,-0.5*Lz,0.5*Lz);CHKERRQ(ierr);
   ierr = IGAAxisInitUniform(axis2,N[2],0.0,Lz,C[2]);CHKERRQ(ierr);
 
-  IGABoundary bnd;
   PetscInt dir=1,side,field;
   for (side=0;side<2;side++) {
-    ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
     for (field=0;field<3;field++) {
-      ierr = IGABoundarySetValue(bnd,field,0.0);CHKERRQ(ierr);
+      ierr = IGASetBoundaryValue(iga,dir,side,field,0.0);CHKERRQ(ierr);
     }
   }
 
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
 
   TS ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);

File demo/Neumann.c

   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   for (dir=0; dir<dim; dir++) {
     for (side=0; side<2; side++) {
-      IGABoundary bnd;
       PetscScalar load = Flux(dir,side);
-      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetLoad(bnd,0,load);CHKERRQ(ierr);
+      ierr = IGASetBoundaryLoad(iga,dir,side,0,load);CHKERRQ(ierr);
     }
   }
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
   ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
   if (!iga->collocation) {
-    ierr = IGASetUserSystem(iga,SystemGalerkin,NULL);CHKERRQ(ierr);
+    ierr = IGASetFormSystem(iga,SystemGalerkin,NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
     ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
   } else {
-    ierr = IGASetUserSystem(iga,SystemCollocation,NULL);CHKERRQ(ierr);
+    ierr = IGASetFormSystem(iga,SystemCollocation,NULL);CHKERRQ(ierr);
   }
 
   {
   }
 
   PetscScalar error;
-  ierr = IGAFormScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(iga,x,1,&error,Error,NULL);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
 
   if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"Error=%G\n",error);CHKERRQ(ierr);}

File demo/PatternFormation.c

   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  ierr = IGASetUserIEFunction(iga,Function,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIEJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIEFunction(iga,Function,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIEJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
 
   PetscReal h  = PetscMin(2.0/N[0],2.0/N[1]);
   PetscReal dt = h/user.delta/15.0;

File demo/PhaseFieldCrystal2D.c

   AppCtx *user = (AppCtx *)mctx;
 
   PetscScalar scalar = 0.0;
-  ierr = IGAFormScalar(user->iga,U,1,&scalar,Stats,mctx);CHKERRQ(ierr);
+  ierr = IGAComputeScalar(user->iga,U,1,&scalar,Stats,mctx);CHKERRQ(ierr);
   PetscReal stats = PetscRealPart(scalar);
 
   PetscReal dt;
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
 
   TS ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);

File demo/Poisson.c

   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   for (dir=0; dir<dim; dir++) {
     for (side=0; side<2; side++) {
-      IGABoundary bnd;
       PetscScalar value = 1.0;
-      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
+      ierr = IGASetBoundaryValue(iga,dir,side,0,value);CHKERRQ(ierr);
     }
   }
 
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
   if (!iga->collocation) {
-    ierr = IGASetUserSystem(iga,SystemGalerkin,NULL);CHKERRQ(ierr);
+    ierr = IGASetFormSystem(iga,SystemGalerkin,NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
   } else {
-    ierr = IGASetUserSystem(iga,SystemCollocation,NULL);CHKERRQ(ierr);
+    ierr = IGASetFormSystem(iga,SystemCollocation,NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_FALSE);CHKERRQ(ierr);
   }
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);

File demo/Poisson1D.c

   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  IGABoundary bnd;
   PetscInt dir=0,side;
   PetscScalar value = 1.0;
   for (side=0; side<2; side++) {
-    ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-    ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);    
+    ierr = IGASetBoundaryValue(iga,dir,side,0,value);CHKERRQ(ierr);
   }
 
   Mat A;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,System,NULL);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 = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
   ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
-  
+
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);

File demo/Poisson2D.c

   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   
-  IGABoundary bnd;
   PetscInt dir,side;
   PetscScalar value = 1.0;
   for (dir=0; dir<2; dir++) {
     for (side=0; side<2; side++) {
-      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
+      ierr = IGASetBoundaryValue(iga,dir,side,0,value);CHKERRQ(ierr);
     }
   }
 
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,System,NULL);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
   
   KSP ksp;

File demo/Poisson3D.c

   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  IGABoundary bnd;
   PetscInt dir,side;
   PetscScalar value = 1.0;
   for (dir=0; dir<3; dir++) {
     for (side=0; side<2; side++) {
-      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
+      ierr = IGASetBoundaryValue(iga,dir,side,0,value);CHKERRQ(ierr);
     }
   }
 
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,System,NULL);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
   
   KSP ksp;

File demo/Richards.c

 #define __FUNCT__ "L2Projection"
 PetscErrorCode L2Projection(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
+  if (p->atboundary) return 0;
+
   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt nen = p->nen;
   Vec b;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,L2Projection,user);CHKERRQ(ierr);
+  ierr = IGASetFormSystem(iga,L2Projection,user);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "ZeroFluxResidual"
+PetscErrorCode ZeroFluxResidual(IGAPoint p,PetscReal dt,
+                                PetscReal shift,const PetscScalar *V,
+                                PetscReal t,const PetscScalar *U,
+                                PetscScalar *R,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+
+  PetscScalar S1[dim];
+  IGAPointFormGrad(p,U,&S1[0]);
+
+  const PetscReal *N0 = (typeof(N0)) p->shape[0];
+  const PetscReal *n = p->normal;
+
+  PetscInt a,i;
+  for (a=0; a<nen; a++) {
+    PetscScalar n_gS = 0.0;
+    for(i=0;i<dim;i++) n_gS += n[i]*S1[i];
+    R[a] = user->penalty*N0[a]*n_gS;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "Residual"
 PetscErrorCode Residual(IGAPoint p,PetscReal dt,
-                           PetscReal shift,const PetscScalar *V,
-                           PetscReal t,const PetscScalar *U,
-                           PetscScalar *R,void *ctx)
+                        PetscReal shift,const PetscScalar *V,
+                        PetscReal t,const PetscScalar *U,
+                        PetscScalar *R,void *ctx)
 {
+  if (p->atboundary)
+    return ZeroFluxResidual(p,dt,shift,V,t,U,R,ctx);
+
   AppCtx *user = (AppCtx *)ctx;
 
   PetscInt a,i;
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "ZeroFluxJacobian"
+PetscErrorCode ZeroFluxJacobian(IGAPoint p,PetscReal dt,
+                                PetscReal shift,const PetscScalar *V,
+                                PetscReal t,const PetscScalar *U,
+                                PetscScalar *J,void *ctx)
+{
+  return 0;
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "Jacobian"
 PetscErrorCode Jacobian(IGAPoint p,PetscReal dt,
                         PetscReal shift,const PetscScalar *V,
                         PetscReal t,const PetscScalar *U,
                         PetscScalar *J,void *ctx)
 {
-
+  if (p->atboundary)
+    return ZeroFluxJacobian(p,dt,shift,V,t,U,J,ctx);
   return 0;
 }
 
-#undef  __FUNCT__
-#define __FUNCT__ "ZeroFluxResidual"
-PetscErrorCode ZeroFluxResidual(IGAPoint p,PetscReal dt,
-                                PetscReal shift,const PetscScalar *V,
-                                PetscReal t,const PetscScalar *U,
-                                PetscScalar *R,void *ctx)
-{
-  AppCtx *user = (AppCtx *)ctx;
-
-  PetscInt nen = p->nen;
-  PetscInt dim = p->dim;
-
-  PetscScalar S1[dim];
-  IGAPointFormGrad(p,U,&S1[0]);
-
-  const PetscReal *N0 = (typeof(N0)) p->shape[0];
-  const PetscReal *n = p->normal;
-
-  PetscInt a,i;
-  for (a=0; a<nen; a++) {
-    PetscScalar n_gS = 0.0;
-    for(i=0;i<dim;i++) n_gS += n[i]*S1[i];
-    R[a] = user->penalty*N0[a]*n_gS;
-  }
-  return 0;
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "ZeroFluxJacobian"
-PetscErrorCode ZeroFluxJacobian(IGAPoint p,PetscReal dt,
-                                PetscReal shift,const PetscScalar *V,
-                                PetscReal t,const PetscScalar *U,
-                                PetscScalar *J,void *ctx)
-{
-
-  return 0;
-}
 
 #undef __FUNCT__
 #define __FUNCT__ "main"
   user.phase_field = PETSC_TRUE;
 
   PetscInt dim = 2, p = 2, k = 1, N = 256, L = 2;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"r_","Richard's Equation Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Richard's Equation Options","IGA");CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   IGA         iga;
   IGAAxis     axis;
-  IGABoundary bnd;
   PetscInt    dir,side;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
   // Boundary conditions
   for(dir=0;dir<dim;dir++)
     for(side=0;side<2;side++){
-      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
       if(dir == dim-1){
-        ierr = IGABoundarySetUserIFunction(bnd,ZeroFluxResidual,&user);CHKERRQ(ierr);
-        ierr = IGABoundarySetUserIJacobian(bnd,ZeroFluxJacobian,&user);CHKERRQ(ierr);
+        ierr = IGASetBoundaryForm(iga,dir,side,PETSC_TRUE);CHKERRQ(ierr);
+        ierr = IGASetBoundaryForm(iga,dir,side,PETSC_TRUE);CHKERRQ(ierr);
         if(side == 0){
-          ierr = IGABoundarySetValue(bnd,0,user.Sin);CHKERRQ(ierr);
+          ierr = IGASetBoundaryValue(iga,dir,side,0,user.Sin);CHKERRQ(ierr);
         }
       }
     }
 
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
 
   TS     ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);

File demo/ShallowWater.c

 
   PetscInt i,j;
   IGAAxis axis;
-  IGABoundary bnd;
   for (i=0; i<2; i++) {
     ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
     ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
     ierr = IGAAxisInitUniform(axis,N[i],0.0,1.0,C[i]);CHKERRQ(ierr);
     if (W[i]) continue;
     for (j=0; j<2; j++) {
-      IGAGetBoundary(iga,i,j,&bnd);
-      IGABoundarySetValue(bnd,1,0.0);
-      IGABoundarySetValue(bnd,2,0.0);
+      IGASetBoundaryValue(iga,i,j,1,0.0);
+      IGASetBoundaryValue(iga,i,j,2,0.0);
     }
   }
 
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
 
   TS ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);

File demo/TwoPhaseTwoComponent.c

 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "LeftInjectionResidual"
+PetscErrorCode LeftInjectionResidual(IGAPoint p,PetscReal dt,
+				     PetscReal shift,const PetscScalar *V,
+				     PetscReal t,const PetscScalar *U,
+				     PetscScalar *R,void *ctx)
+{
+  PetscInt a,nen;
+  IGAPointGetSizes(p,&nen,0,0);
+
+  const PetscReal *N0;
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+
+  PetscScalar Qh = 0;
+  if(t <= 5e5) Qh = -5.57e-6; // inflow
+
+  PetscScalar (*Re)[2] = (PetscScalar (*)[2])R;
+  for (a=0; a<nen; a++) {
+    Re[a][0] = 0.0;
+    Re[a][1] = N0[a]*Qh;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "Residual"
 PetscErrorCode Residual(IGAPoint p,PetscReal dt,
 			PetscReal shift,const PetscScalar *V,
 			PetscReal t,const PetscScalar *U,
 			PetscScalar *R,void *ctx)
 {
+  if (p->atboundary)
+    return LeftInjectionResidual(p,dt,shift,V,t,U,R,ctx);
+
   AppCtx *user = (AppCtx *)ctx;
   PetscScalar rholw = user->rholw; 
   PetscScalar porosity = user->porosity; 
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "Jacobian"
-PetscErrorCode Jacobian(IGAPoint p,PetscReal dt,
-			PetscReal shift,const PetscScalar *V,
-			PetscReal t,const PetscScalar *U,
-			PetscScalar *J,void *ctx)
-{
-  // for now use the option -snes_fd_color for the Jacobian
-  return 0;
-}
-
-#undef  __FUNCT__
 #define __FUNCT__ "LeftInjectionJacobian"
 PetscErrorCode LeftInjectionJacobian(IGAPoint p,PetscReal dt,
 				     PetscReal shift,const PetscScalar *V,
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "LeftInjectionResidual"
-PetscErrorCode LeftInjectionResidual(IGAPoint p,PetscReal dt,
-				     PetscReal shift,const PetscScalar *V,
-				     PetscReal t,const PetscScalar *U,
-				     PetscScalar *R,void *ctx)
+#define __FUNCT__ "Jacobian"
+PetscErrorCode Jacobian(IGAPoint p,PetscReal dt,
+			PetscReal shift,const PetscScalar *V,
+			PetscReal t,const PetscScalar *U,
+			PetscScalar *J,void *ctx)
 {
-  PetscInt a,nen;
-  IGAPointGetSizes(p,&nen,0,0);
-
-  const PetscReal *N0;
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
-
-  PetscScalar Qh = 0;
-  if(t <= 5e5) Qh = -5.57e-6; // inflow
-
-  PetscScalar (*Re)[2] = (PetscScalar (*)[2])R;
-  for (a=0; a<nen; a++) {
-    Re[a][0] = 0.0;
-    Re[a][1] = N0[a]*Qh;
-  }
+  if (p->atboundary)
+    return LeftInjectionJacobian(p,dt,shift,V,t,U,J,ctx);
+  // for now use the option -snes_fd_color for the Jacobian
   return 0;
 }
 
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   user.iga = iga;
 
-  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetFormIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
 
-  IGABoundary bnd;
-  ierr = IGAGetBoundary(iga,0,0,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetUserIFunction(bnd,LeftInjectionResidual,&user);CHKERRQ(ierr);
-  ierr = IGABoundarySetUserIJacobian(bnd,LeftInjectionJacobian,&user);CHKERRQ(ierr);
-  ierr = IGAGetBoundary(iga,0,1,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,10.0);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,1, 0.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryForm(iga,0,0,PETSC_TRUE);CHKERRQ(ierr);
+  ierr = IGASetBoundaryForm(iga,0,0,PETSC_TRUE);CHKERRQ(ierr);
+
+  ierr = IGASetBoundaryValue(iga,0,1,0,10.0);CHKERRQ(ierr);
+  ierr = IGASetBoundaryValue(iga,0,1,1, 0.0);CHKERRQ(ierr);
 
   TS     ts;
   ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);

File include/petiga.h

 
 /* ---------------------------------------------------------------- */
 
+typedef struct _p_IGA         *IGA;
 typedef struct _n_IGAAxis     *IGAAxis;
-typedef struct _n_IGABoundary *IGABoundary;
 typedef struct _n_IGARule     *IGARule;
 typedef struct _n_IGABasis    *IGABasis;
+typedef struct _n_IGAForm     *IGAForm;
 typedef struct _n_IGAElement  *IGAElement;
 typedef struct _n_IGAPoint    *IGAPoint;
 
   PetscReal  bnd_point[2];
   PetscReal *bnd_value[2]; /* [nen][d+1] */
 };
+
 PETSC_EXTERN PetscErrorCode IGABasisCreate(IGABasis *basis);
 PETSC_EXTERN PetscErrorCode IGABasisDestroy(IGABasis *basis);
 PETSC_EXTERN PetscErrorCode IGABasisReset(IGABasis basis);
 
 /* ---------------------------------------------------------------- */
 
-typedef PetscErrorCode (*IGAUserScalar)    (IGAPoint point,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx);
-typedef PetscErrorCode (*IGAUserVector)    (IGAPoint point,PetscScalar *F,void *ctx);
-typedef PetscErrorCode (*IGAUserMatrix)    (IGAPoint point,PetscScalar *K,void *ctx);
-typedef PetscErrorCode (*IGAUserSystem)    (IGAPoint point,PetscScalar *K,PetscScalar *F,void *ctx);
-typedef PetscErrorCode (*IGAUserFunction)  (IGAPoint point,const PetscScalar *U,PetscScalar *F,void *ctx);
-typedef PetscErrorCode (*IGAUserJacobian)  (IGAPoint point,const PetscScalar *U,PetscScalar *J,void *ctx);
-typedef PetscErrorCode (*IGAUserIFunction) (IGAPoint point,PetscReal dt,
+typedef PetscErrorCode (*IGAFormScalar)    (IGAPoint point,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx);
+typedef PetscErrorCode (*IGAFormVector)    (IGAPoint point,PetscScalar *F,void *ctx);
+typedef PetscErrorCode (*IGAFormMatrix)    (IGAPoint point,PetscScalar *K,void *ctx);
+typedef PetscErrorCode (*IGAFormSystem)    (IGAPoint point,PetscScalar *K,PetscScalar *F,void *ctx);
+typedef PetscErrorCode (*IGAFormFunction)  (IGAPoint point,const PetscScalar *U,PetscScalar *F,void *ctx);
+typedef PetscErrorCode (*IGAFormJacobian)  (IGAPoint point,const PetscScalar *U,PetscScalar *J,void *ctx);
+typedef PetscErrorCode (*IGAFormIFunction) (IGAPoint point,PetscReal dt,
                                             PetscReal a,const PetscScalar *V,
                                             PetscReal t,const PetscScalar *U,
                                             PetscScalar *F,void *ctx);
-typedef PetscErrorCode (*IGAUserIJacobian) (IGAPoint point,PetscReal dt,
+typedef PetscErrorCode (*IGAFormIJacobian) (IGAPoint point,PetscReal dt,
                                             PetscReal a,const PetscScalar *V,
                                             PetscReal t,const PetscScalar *U,
                                             PetscScalar *J,void *ctx);
-typedef PetscErrorCode (*IGAUserIFunction2)(IGAPoint point,PetscReal dt,
+typedef PetscErrorCode (*IGAFormIFunction2)(IGAPoint point,PetscReal dt,
                                             PetscReal a,const PetscScalar *A,
                                             PetscReal v,const PetscScalar *V,
                                             PetscReal t,const PetscScalar *U,
                                             PetscScalar *F,void *ctx);
-typedef PetscErrorCode (*IGAUserIJacobian2)(IGAPoint point,PetscReal dt,
+typedef PetscErrorCode (*IGAFormIJacobian2)(IGAPoint point,PetscReal dt,
                                             PetscReal a,const PetscScalar *A,
                                             PetscReal v,const PetscScalar *V,
                                             PetscReal t,const PetscScalar *U,
                                             PetscScalar *J,void *ctx);
-typedef PetscErrorCode (*IGAUserIEFunction)(IGAPoint point,PetscReal dt,
+typedef PetscErrorCode (*IGAFormIEFunction)(IGAPoint point,PetscReal dt,
                                             PetscReal a,const PetscScalar *V,
                                             PetscReal t,const PetscScalar *U,
                                             PetscReal t0,const PetscScalar *U0,
                                             PetscScalar *F,void *ctx);
-typedef PetscErrorCode (*IGAUserIEJacobian)(IGAPoint point,PetscReal dt,
+typedef PetscErrorCode (*IGAFormIEJacobian)(IGAPoint point,PetscReal dt,
                                             PetscReal a,const PetscScalar *V,
                                             PetscReal t,const PetscScalar *U,
                                             PetscReal t0,const PetscScalar *U0,
                                             PetscScalar *J,void *ctx);
 
-typedef struct _IGAUserOps *IGAUserOps;
-struct _IGAUserOps {
+typedef struct _IGAFormBC *IGAFormBC;
+struct _IGAFormBC {
+  PetscInt    count;
+  PetscInt    field[64];
+  PetscScalar value[64];
+};
+
+typedef struct _IGAFormOps *IGAFormOps;
+struct _IGAFormOps {
   /**/
-  IGAUserVector     Vector;
+  IGAFormVector     Vector;
   void              *VecCtx;
-  IGAUserMatrix     Matrix;
+  IGAFormMatrix     Matrix;
   void              *MatCtx;
-  IGAUserSystem     System;
+  IGAFormSystem     System;
   void              *SysCtx;
   /**/
-  IGAUserFunction   Function;
+  IGAFormFunction   Function;
   void              *FunCtx;
-  IGAUserJacobian   Jacobian;
+  IGAFormJacobian   Jacobian;
   void              *JacCtx;
   /**/
-  IGAUserIFunction  IFunction;
-  IGAUserIFunction2 IFunction2;
+  IGAFormIFunction  IFunction;
+  IGAFormIFunction2 IFunction2;
   void              *IFunCtx;
-  IGAUserIJacobian  IJacobian;
-  IGAUserIJacobian2 IJacobian2;
+  IGAFormIJacobian  IJacobian;
+  IGAFormIJacobian2 IJacobian2;
   void              *IJacCtx;
   /**/
-  IGAUserIEFunction IEFunction;
+  IGAFormIEFunction IEFunction;
   void              *IEFunCtx;
-  IGAUserIEJacobian IEJacobian;
+  IGAFormIEJacobian IEJacobian;
   void              *IEJacCtx;
 };
 
-/* ---------------------------------------------------------------- */
-
-struct _n_IGABoundary {
+struct _n_IGAForm {
   PetscInt refct;
   /**/
-  PetscInt    dof;
-  /**/
-  PetscInt    count;
-  PetscInt    *field;
-  PetscScalar *value;
-  /**/
-  PetscInt    nload;
-  PetscInt    *iload;
-  PetscScalar *vload;
-  /**/
-  IGAUserOps userops;
+  IGAFormOps ops;
+  PetscInt   dof;
+  IGAFormBC  value[3][2];
+  IGAFormBC  load [3][2];
+  PetscBool  visit[3][2];
 };
-PETSC_EXTERN PetscErrorCode IGABoundaryCreate(IGABoundary *boundary);
-PETSC_EXTERN PetscErrorCode IGABoundaryDestroy(IGABoundary *boundary);
-PETSC_EXTERN PetscErrorCode IGABoundaryReset(IGABoundary boundary);
-PETSC_EXTERN PetscErrorCode IGABoundaryReference(IGABoundary boundary);
-PETSC_EXTERN PetscErrorCode IGABoundaryInit(IGABoundary boundary,PetscInt dof);
-PETSC_EXTERN PetscErrorCode IGABoundaryClear(IGABoundary boundary);
-PETSC_EXTERN PetscErrorCode IGABoundarySetValue(IGABoundary boundary,PetscInt field,PetscScalar value);
-PETSC_EXTERN PetscErrorCode IGABoundarySetLoad (IGABoundary boundary,PetscInt field,PetscScalar value);
 
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserVector    (IGABoundary boundary,IGAUserVector     Vector,    void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserMatrix    (IGABoundary boundary,IGAUserMatrix     Matrix,    void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserSystem    (IGABoundary boundary,IGAUserSystem     System,    void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserFunction  (IGABoundary boundary,IGAUserFunction   Function,  void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserJacobian  (IGABoundary boundary,IGAUserJacobian   Jacobian,  void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserIFunction (IGABoundary boundary,IGAUserIFunction  IFunction, void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserIJacobian (IGABoundary boundary,IGAUserIJacobian  IJacobian, void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserIFunction2(IGABoundary boundary,IGAUserIFunction2 IFunction, void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserIJacobian2(IGABoundary boundary,IGAUserIJacobian2 IJacobian, void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserIEFunction(IGABoundary boundary,IGAUserIEFunction IEFunction,void *ctx);
-PETSC_EXTERN PetscErrorCode IGABoundarySetUserIEJacobian(IGABoundary boundary,IGAUserIEJacobian IEJacobian,void *ctx);
+PETSC_EXTERN PetscErrorCode IGAGetForm(IGA iga,IGAForm *form);
+PETSC_EXTERN PetscErrorCode IGASetForm(IGA iga,IGAForm form);
 
+PETSC_EXTERN PetscErrorCode IGAFormCreate(IGAForm *form);
+PETSC_EXTERN PetscErrorCode IGAFormDestroy(IGAForm *form);
+PETSC_EXTERN PetscErrorCode IGAFormReset(IGAForm form);
+PETSC_EXTERN PetscErrorCode IGAFormReference(IGAForm form);
 
-typedef struct _p_IGA *IGA;
+PETSC_EXTERN PetscErrorCode IGAFormSetBoundaryValue(IGAForm form,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
+PETSC_EXTERN PetscErrorCode IGAFormSetBoundaryLoad (IGAForm form,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
+PETSC_EXTERN PetscErrorCode IGAFormSetBoundaryForm (IGAForm form,PetscInt axis,PetscInt side,PetscBool flag);
+PETSC_EXTERN PetscErrorCode IGAFormClearBoundary   (IGAForm form,PetscInt axis,PetscInt side);
+
+PETSC_EXTERN PetscErrorCode IGAFormSetVector    (IGAForm form,IGAFormVector     Vector,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetMatrix    (IGAForm form,IGAFormMatrix     Matrix,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetSystem    (IGAForm form,IGAFormSystem     System,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetFunction  (IGAForm form,IGAFormFunction   Function,  void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetJacobian  (IGAForm form,IGAFormJacobian   Jacobian,  void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetIFunction (IGAForm form,IGAFormIFunction  IFunction, void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetIJacobian (IGAForm form,IGAFormIJacobian  IJacobian, void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetIFunction2(IGAForm form,IGAFormIFunction2 IFunction, void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetIJacobian2(IGAForm form,IGAFormIJacobian2 IJacobian, void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetIEFunction(IGAForm form,IGAFormIEFunction IEFunction,void *ctx);
+PETSC_EXTERN PetscErrorCode IGAFormSetIEJacobian(IGAForm form,IGAFormIEJacobian IEJacobian,void *ctx);
+
+PETSC_EXTERN PetscErrorCode IGASetBoundaryValue(IGA iga,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
+PETSC_EXTERN PetscErrorCode IGASetBoundaryLoad (IGA iga,PetscInt axis,PetscInt side,PetscInt field,PetscScalar value);
+PETSC_EXTERN PetscErrorCode IGASetBoundaryForm (IGA iga,PetscInt axis,PetscInt side,PetscBool flag);
+
+PETSC_EXTERN PetscErrorCode IGASetFormVector    (IGA iga,IGAFormVector     Vector,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormMatrix    (IGA iga,IGAFormMatrix     Matrix,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormSystem    (IGA iga,IGAFormSystem     System,    void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormFunction  (IGA iga,IGAFormFunction   Function,  void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormJacobian  (IGA iga,IGAFormJacobian   Jacobian,  void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormIFunction (IGA iga,IGAFormIFunction  IFunction, void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormIJacobian (IGA iga,IGAFormIJacobian  IJacobian, void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormIFunction2(IGA iga,IGAFormIFunction2 IFunction, void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormIJacobian2(IGA iga,IGAFormIJacobian2 IJacobian, void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormIEFunction(IGA iga,IGAFormIEFunction IEFunction,void *ctx);
+PETSC_EXTERN PetscErrorCode IGASetFormIEJacobian(IGA iga,IGAFormIEJacobian IEJacobian,void *ctx);
+
+/* ---------------------------------------------------------------- */
 
 typedef struct _IGAOps *IGAOps;
 struct _IGAOps {
   PetscInt   setupstage;
   PetscBool  collocation;
 
-  IGAUserOps userops;
   VecType    vectype;
   MatType    mattype;
   char       **fieldname;
 
   IGAAxis     axis[3];
   IGARule     rule[3];
-  IGABoundary boundary[3][2];
   IGABasis    basis[3];
+  IGAForm     form;
 
   IGAElement  iterator;
 
 PETSC_EXTERN PetscErrorCode IGAGetAxis(IGA iga,PetscInt i,IGAAxis *axis);
 PETSC_EXTERN PetscErrorCode IGAGetRule(IGA iga,PetscInt i,IGARule *rule);
 PETSC_EXTERN PetscErrorCode IGAGetBasis(IGA iga,PetscInt i,IGABasis *basis);
-PETSC_EXTERN PetscErrorCode IGAGetBoundary(IGA iga,PetscInt i,PetscInt side,IGABoundary *boundary);
 
 PETSC_EXTERN PetscErrorCode IGAGetComm(IGA iga,MPI_Comm *comm);
 
 PETSC_EXTERN PetscErrorCode IGAGetLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);
 PETSC_EXTERN PetscErrorCode IGARestoreLocalVecArray(IGA iga,Vec gvec,Vec *lvec,const PetscScalar *array[]);
 
-PETSC_EXTERN PetscErrorCode IGASetUserVector    (IGA iga,IGAUserVector     Vector,    void *ctx);
-PETSC_EXTERN PetscErrorCode IGASetUserMatrix    (IGA iga,IGAUserMatrix     Matrix,    void *ctx);
-PETSC_EXTERN PetscErrorCode IGASetUserSystem    (IGA iga,IGAUserSystem     System,    void *ctx);
-PETSC_EXTERN PetscErrorCode IGASetUserFunction  (IGA iga,IGAUserFunction   Function,  void *ctx);
-PETSC_EXTERN PetscErrorCode IGASetUserJacobian  (IGA iga,IGAUserJacobian   Jacobian,  void *ctx);
-PETSC_EXTERN PetscErrorCode IGASetUserIFunction (IGA iga,IGAUserIFunction  IFunction, void *ctx);
-PETSC_EXTERN PetscErrorCode IGASetUserIJacobian (IGA iga,IGAUserIJacobian  IJacobian, void *ctx);
-PETSC_EXTERN PetscErrorCode IGASetUserIFunction2(IGA iga,IGAUserIFunction2 IFunction, void *ctx);
-PETSC_EXTERN PetscErrorCode IGASetUserIJacobian2(IGA iga,IGAUserIJacobian2 IJacobian, void *ctx);
-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 PetscErrorCode DMIGASetIGA(DM dm,IGA iga);
 PETSC_EXTERN PetscErrorCode DMIGAGetIGA(DM dm,IGA *iga);
 
+
+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[]);
+
 /* ---------------------------------------------------------------- */
 
 struct _n_IGAElement {
   PetscInt  width[3];
   PetscInt  sizes[3];
   PetscInt  ID[3];
-  IGABasis  BD[3];
+  /**/
   PetscBool atboundary;
   PetscInt  boundary_id;
   /**/
 PETSC_EXTERN PetscErrorCode IGABeginElement(IGA iga,IGAElement *element);
 PETSC_EXTERN PetscBool      IGANextElement(IGA iga,IGAElement element);
 PETSC_EXTERN PetscErrorCode IGAEndElement(IGA iga,IGAElement *element);
-PETSC_EXTERN PetscBool      IGAElementNextUserOps(IGAElement element,IGAUserOps *userops);
+PETSC_EXTERN PetscBool      IGAElementNextForm(IGAElement element,PetscBool visit[][2]);
 PETSC_EXTERN PetscErrorCode IGAElementGetPoint(IGAElement element,IGAPoint *point);
 PETSC_EXTERN PetscErrorCode IGAElementBeginPoint(IGAElement element,IGAPoint *point);
 PETSC_EXTERN PetscBool      IGAElementNextPoint(IGAElement element,IGAPoint point);
 PETSC_EXTERN PetscErrorCode IGAPointGetParent(IGAPoint point,IGAElement *element);
 PETSC_EXTERN PetscErrorCode IGAPointGetIndex(IGAPoint point,PetscInt *index);
 PETSC_EXTERN PetscErrorCode IGAPointGetCount(IGAPoint point,PetscInt *count);
+PETSC_EXTERN PetscErrorCode IGAPointAtBoundary(IGAPoint point,PetscInt *axis,PetscInt *side);
 PETSC_EXTERN PetscErrorCode IGAPointGetSizes(IGAPoint point,PetscInt *neq,PetscInt *nen,PetscInt *dof);
 PETSC_EXTERN PetscErrorCode IGAPointGetDims(IGAPoint point,PetscInt *dim,PetscInt *nsd,PetscInt *npd);
 PETSC_EXTERN PetscErrorCode IGAPointGetQuadrature(IGAPoint point,PetscReal *weigth,PetscReal *detJac);
 PETSC_EXTERN PetscLogEvent IGA_FormIFunction;
 PETSC_EXTERN PetscLogEvent IGA_FormIJacobian;
 
-PETSC_EXTERN PetscErrorCode IGAFormScalar(IGA iga,Vec U,PetscInt n,PetscScalar S[],
-                                          IGAUserScalar Scalar,void *ctx);
+PETSC_EXTERN PetscErrorCode IGAComputeScalar(IGA iga,Vec U,
+                                             PetscInt n,PetscScalar S[],
+                                             IGAFormScalar Scalar,void *ctx);
 
 #define PCIGAEBE "igaebe"
 #define PCIGABBB "igabbb"
 #endif
 
 #if defined(PETSC_USE_DEBUG)
-#define IGACheckUserOp(iga,arg,UserOp) do {             \
-    if (!iga->userops->UserOp)                          \
+#define IGACheckFormOp(iga,arg,FormOp) do {             \
+    if (!iga->form->ops->FormOp)                        \
       SETERRQ4(((PetscObject)iga)->comm,PETSC_ERR_USER, \
-               "Must call IGASetUser%s() "              \
+               "Must call IGASetForm%s() "              \
                "on argument %D \"%s\" before %s()",     \
-               #UserOp,(arg),#iga,PETSC_FUNCTION_NAME); \
+               #FormOp,(arg),#iga,PETSC_FUNCTION_NAME); \
     } while (0)
 #else
-#define IGACheckUserOp(iga,arg,UserOp) do {} while (0)
+#define IGACheckFormOp(iga,arg,FormOp) do {} while (0)
 #endif
 
 /* ---------------------------------------------------------------- */

File src/makefile

 petigaaxis.c \
 petigarule.c \
 petigabasis.c \
-petigabound.c \
+petigaform.c \
 petigaelem.c \
 petigapoint.c \
 petigavec.c \
 petigamat.c \
 petigadm.c \
-petigascl.c \
+petigacomp.c \
 petigapcb.c \
 petigapce.c \
 petigapc.c \

File src/petiga.c

 
   *_iga = iga;
 
-  ierr = PetscNew(struct _IGAUserOps,&iga->userops);CHKERRQ(ierr);
   iga->vectype = NULL;
   iga->mattype = NULL;
 
   iga->dof = -1;
   iga->order = -1;
 
+  for (i=0; i<3; i++)
+    iga->proc_sizes[i] = -1;
+
   for (i=0; i<3; i++) {
     ierr = IGAAxisCreate(&iga->axis[i]);CHKERRQ(ierr);
     ierr = IGARuleCreate(&iga->rule[i]);CHKERRQ(ierr);
     ierr = IGABasisCreate(&iga->basis[i]);CHKERRQ(ierr);
-    ierr = IGABoundaryCreate(&iga->boundary[i][0]);CHKERRQ(ierr);
-    ierr = IGABoundaryCreate(&iga->boundary[i][1]);CHKERRQ(ierr);
-    iga->proc_sizes[i] = -1;
   }
+  ierr = IGAFormCreate(&iga->form);CHKERRQ(ierr);
   ierr = IGAElementCreate(&iga->iterator);CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   if (--((PetscObject)iga)->refct > 0) PetscFunctionReturn(0);
 
-  ierr = PetscFree(iga->userops);CHKERRQ(ierr);
   ierr = PetscFree(iga->vectype);CHKERRQ(ierr);
   ierr = PetscFree(iga->mattype);CHKERRQ(ierr);
   if (iga->fieldname) {
     ierr = IGAAxisDestroy(&iga->axis[i]);CHKERRQ(ierr);
     ierr = IGARuleDestroy(&iga->rule[i]);CHKERRQ(ierr);
     ierr = IGABasisDestroy(&iga->basis[i]);CHKERRQ(ierr);
-    ierr = IGABoundaryDestroy(&iga->boundary[i][0]);CHKERRQ(ierr);
-    ierr = IGABoundaryDestroy(&iga->boundary[i][1]);CHKERRQ(ierr);
   }
+  ierr = IGAFormDestroy(&iga->form);CHKERRQ(ierr);
   ierr = IGAElementDestroy(&iga->iterator);CHKERRQ(ierr);
 
   ierr = IGAReset(iga);CHKERRQ(ierr);
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "IGAGetBoundary"
-/*@
-   IGAGetBoundary - Returns a pointer to a specific side of the i^th
-   boundary associated with the IGA.
-
-   Not Collective
-
-   Input Parameters:
-+  iga - the IGA context
-.  i - the boundary index
--  side - the side index: 0 (left) or 1 (right)
-
-   Output Parameter:
-.  boundary - the boundary context
-
-   Notes:
-   A side marker of 0 corresponds to the boundary associated to the
-   minimum knot value of the i^th axis. A side marker of 1 corresponds
-   to the boundary associated to the maximum knot value of the i^th
-   axis.
-
-   Level: normal
-
-.keywords: IGA, boundary
-@*/
-PetscErrorCode IGAGetBoundary(IGA iga,PetscInt i,PetscInt side,IGABoundary *boundary)
-{
-  PetscInt       dim;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  PetscValidPointer(boundary,4);
-  dim = (iga->dim > 0) ? iga->dim : 3;
-  if (i < 0)    SETERRQ1(((PetscObject)iga)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Index %D must be nonnegative",i);
-  if (i >= dim) SETERRQ2(((PetscObject)iga)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Index %D, but dimension %D",i,iga->dim);
-  if (iga->dof <= 0) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call IGASetDof() first");
-  if (side < 0) side = 0; /* XXX error ?*/
-  if (side > 1) side = 1; /* XXX error ?*/
-  if (iga->boundary[i][side]->dof != iga->dof) {
-    ierr = IGABoundaryInit(iga->boundary[i][side],iga->dof);CHKERRQ(ierr);
-  }
-  *boundary = iga->boundary[i][side];
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
 #define __FUNCT__ "IGAGetBasis"
 PetscErrorCode IGAGetBasis(IGA iga,PetscInt i,IGABasis *basis)
 {
   /* --- Stage 3 --- */
   iga->setupstage = 3;
 
-  for (i=iga->dim; i<3; i++) {
-    ierr = IGABoundaryReset(iga->boundary[i][0]);CHKERRQ(ierr);
-    ierr = IGABoundaryReset(iga->boundary[i][1]);CHKERRQ(ierr);
-  }
-
   if (iga->order < 0)
     for (i=0; i<iga->dim; i++)
       iga->order = PetscMax(iga->order,iga->axis[i]->p);
   iga->mattype = mtype;
   PetscFunctionReturn(0);
 }
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserVector"
-/*@
-   IGASetUserSystem - Set the user callback to form the vector
-   which represents the discretized L(w).
-
-   Logically collective on IGA
-
-   Input Parameters:
-+  iga - the IGA context
-.  Vector - the function which evaluates L(w)
--  ctx - user-defined context for evaluation routine (may be NULL)
-
-   Details of Vector:
-$  PetscErrorCode Vector(IGAPoint p,PetscScalar *F,void *ctx);
-
-+  p - point at which to evaluate L(w)
-.  F - contribution to L(w)
--  ctx - user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, setup linear system, vector assembly
-@*/
-PetscErrorCode IGASetUserVector(IGA iga,IGAUserVector Vector,void *VecCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (Vector) iga->userops->Vector = Vector;
-  if (VecCtx) iga->userops->VecCtx = VecCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserMatrix"
-/*@
-   IGASetUserSystem - Set the user callback to form the matrix and vector
-   which represents the discretized a(w,u).
-
-   Logically collective on IGA
-
-   Input Parameters:
-+  iga - the IGA context
-.  Matrix - the function which evaluates a(w,u)
--  ctx - user-defined context for evaluation routine (may be NULL)
-
-   Details of System:
-$  PetscErrorCode System(IGAPoint p,PetscScalar *K,void *ctx);
-
-+  p - point at which to evaluate a(w,u)
-.  K - contribution to a(w,u)
--  ctx - user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, setup linear system, matrix assembly
-@*/
-PetscErrorCode IGASetUserMatrix(IGA iga,IGAUserMatrix Matrix,void *MatCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (Matrix) iga->userops->Matrix = Matrix;
-  if (MatCtx) iga->userops->MatCtx = MatCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserSystem"
-/*@
-   IGASetUserSystem - Set the user callback to form the matrix and vector
-   which represents the discretized a(w,u) = L(w).
-
-   Logically collective on IGA
-
-   Input Parameters:
-+  iga - the IGA context
-.  System - the function which evaluates a(w,u) and L(w)
--  ctx - user-defined context for evaluation routine (may be NULL)
-
-   Details of System:
-$  PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx);
-
-+  p - point at which to evaluate a(w,u)=L(w)
-.  K - contribution to a(w,u)
-.  F - contribution to L(w)
--  ctx - user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, setup linear system, matrix assembly, vector assembly
-@*/
-PetscErrorCode IGASetUserSystem(IGA iga,IGAUserSystem System,void *SysCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (System) iga->userops->System = System;
-  if (SysCtx) iga->userops->SysCtx = SysCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserFunction"
-PetscErrorCode IGASetUserFunction(IGA iga,IGAUserFunction Function,void *FunCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (Function) iga->userops->Function = Function;
-  if (FunCtx)   iga->userops->FunCtx   = FunCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserJacobian"
-PetscErrorCode IGASetUserJacobian(IGA iga,IGAUserJacobian Jacobian,void *JacCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (Jacobian) iga->userops->Jacobian = Jacobian;
-  if (JacCtx)   iga->userops->JacCtx   = JacCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserIFunction"
-/*@
-   IGASetUserIFunction - Set the function which computes the residual
-   R(u)=0 for use with implicit time stepping routines.
-
-   Logically Collective on TS
-
-   Input Parameter:
-+  iga - the IGA context
-.  IFunction - the function evaluation routine
--  IFunCtx - user-defined context for private data for the function evaluation routine (may be NULL)
-
-   Details of IFunction:
-$  PetscErrorCode IFunction(IGAPoint p,PetscReal dt,
-                            PetscReal shift,const PetscScalar *V,
-                            PetscReal t,const PetscScalar *U,
-                            PetscScalar *R,void *ctx);
-
-+  p - point at which to compute the residual
-.  dt - time step size
-.  shift - positive parameter which depends on the time integration method (XXX Should this be here?)
-.  V - time derivative of the state vector
-.  t - time at step/stage being solved
-.  U - state vector
-.  R - function vector
--  ctx - [optional] user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, options
-@*/
-PetscErrorCode IGASetUserIFunction(IGA iga,IGAUserIFunction IFunction,void *IFunCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (IFunction) iga->userops->IFunction = IFunction;
-  if (IFunCtx)   iga->userops->IFunCtx   = IFunCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserIJacobian"
-/*@
-   IGASetUserIJacobian - Set the function to compute the matrix dF/dU
-   + shift*dF/dU_t where F(t,U,U_t) is the function you provided with
-   IGASetUserIFunction().
-
-   Logically Collective on TS
-
-   Input Parameter:
-+  iga - the IGA context
-.  IJacobian - the Jacobian evaluation routine
--  IJacCtx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
-
-   Details of IJacobian:
-$  PetscErrorCode IJacobian(IGAPoint p,PetscReal dt,
-                            PetscReal shift,const PetscScalar *V,
-                            PetscReal t,const PetscScalar *U,
-                            PetscScalar *J,void *ctx);
-
-+  p - point at which to compute the Jacobian
-.  dt - time step size
-.  shift - positive parameter which depends on the time integration method
-.  V - time derivative of the state vector
-.  t - time at step/stage being solved
-.  U - state vector
-.  J - Jacobian matrix
--  ctx - [optional] user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, options
-@*/
-PetscErrorCode IGASetUserIJacobian(IGA iga,IGAUserIJacobian IJacobian,void *IJacCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (IJacobian) iga->userops->IJacobian = IJacobian;
-  if (IJacCtx)   iga->userops->IJacCtx   = IJacCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserIFunction2"
-/*@
-   IGASetUserIFunction - Set the function which computes the residual
-       F(t,U_tt,U_t,U)=0 for use with implicit time stepping routines.
-
-   Logically Collective on TS
-
-   Input Parameter:
-+  iga - the IGA context
-.  IFunction - the function evaluation routine
--  IFunCtx - user-defined context for private data for the function evaluation routine (may be NULL)
-
-   Details of IFunction:
-$  PetscErrorCode IFunction(IGAPoint p,PetscReal dt,
-                            PetscReal a,const PetscScalar *A,
-                            PetscReal v,const PetscScalar *V,
-                            PetscReal t,const PetscScalar *U,
-                            PetscScalar *F,void *ctx);
-
-+  p - point at which to compute the residual
-.  dt - time step size
-.  a - positive parameter which depends on the time integration method
-.  A - second time derivative of the state vector
-.  v - positive parameter which depends on the time integration method
-.  V - time derivative of the state vector
-.  t - time at step/stage being solved
-.  U - state vector
-.  F - function vector
--  ctx - [optional] user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, options
-@*/
-PetscErrorCode IGASetUserIFunction2(IGA iga,IGAUserIFunction2 IFunction,void *IFunCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (IFunction) iga->userops->IFunction2 = IFunction;
-  if (IFunCtx)   iga->userops->IFunCtx    = IFunCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserIJacobian2"
-/*@
-   IGASetUserIJacobian2 - Set the function to compute the matrix
-       J = a*dF/dU_tt + v*dF/dU_t + dF/dU  where F(t,U_tt,U_t,U) is
-       the function you provided with IGASetUserIFunction2().
-
-   Logically Collective on TS
-
-   Input Parameter:
-+  iga       - the IGA context
-.  IJacobian - the Jacobian evaluation routine
--  IJacCtx   - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
-
-   Details of IJacobian:
-$  PetscErrorCode IJacobian(IGAPoint p,PetscReal dt,
-                            PetscReal a,const PetscScalar *A,
-                            PetscReal v,const PetscScalar *V,
-                            PetscReal t,const PetscScalar *U,
-                            PetscScalar *J,void *ctx);
-
-+  p   - point at which to compute the Jacobian
-.  dt  - time step size
-.  a   - positive parameter which depends on the time integration method
-.  A   - second time derivative of the state vector
-.  v   - positive parameter which depends on the time integration method
-.  V   - time derivative of the state vector
-.  t   - time at step/stage being solved
-.  U   - state vector
-.  J   - Jacobian matrix
--  ctx - [optional] user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, options
-@*/
-PetscErrorCode IGASetUserIJacobian2(IGA iga,IGAUserIJacobian2 IJacobian,void *IJacCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (IJacobian) iga->userops->IJacobian2 = IJacobian;
-  if (IJacCtx)   iga->userops->IJacCtx    = IJacCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserIEFunction"
-/*@
-   IGASetUserIEFunction - Set the function which computes the residual
-   R(u)=0 for use with explicit or implicit time stepping routines.
-
-   Logically Collective on TS
-
-   Input Parameter:
-+  iga - the IGA context
-.  IEFunction - the function evaluation routine
--  IEFunCtx - user-defined context for private data for the function evaluation routine (may be NULL)
-
-   Details of IEFunction:
-$  PetscErrorCode IEFunction(IGAPoint p,PetscReal dt,
-                             PetscReal shift,const PetscScalar *V0,
-                             PetscReal t1,const PetscScalar *U1,
-                             PetscReal t0,const PetscScalar *U0,
-                             PetscScalar *R,void *ctx);
-
-+  p - point at which to compute the residual
-.  dt - time step size
-.  shift - positive parameter which depends on the time integration method (XXX Should this be here?)
-.  V0 - time derivative of the state vector at t0
-.  t1 - time at step/stage being solved
-.  U1 - state vector at t1
-.  t0 - time at current step
-.  U0 - state vector at t0
-.  R - function vector
--  ctx - [optional] user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, options
-@*/
-PetscErrorCode IGASetUserIEFunction(IGA iga,IGAUserIEFunction IEFunction,void *IEFunCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (IEFunction) iga->userops->IEFunction = IEFunction;
-  if (IEFunCtx)   iga->userops->IEFunCtx   = IEFunCtx;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGASetUserIEJacobian"
-/*@
-   IGASetUserIEJacobian - Set the function to compute the matrix dF/dU
-   + shift*dF/dU_t where F(t,U,U_t) is the function you provided with
-   IGASetUserIEFunction(). For use with implicit or explicit TS methods.
-
-   Logically Collective on TS
-
-   Input Parameter:
-+  iga - the IGA context
-.  IEJacobian - the Jacobian evaluation routine
--  IEJacCtx - user-defined context for private data for the Jacobian evaluation routine (may be NULL)
-
-   Details of IEJacobian:
-$  PetscErrorCode IEJacobian(IGAPoint p,PetscReal dt,
-                             PetscReal shift,const PetscScalar *V0,
-                             PetscReal t1,const PetscScalar *U1,
-                             PetscReal t0,const PetscScalar *U0,
-                             PetscScalar *J,void *ctx);
-
-+  p - point at which to compute the Jacobian
-.  dt - time step size
-.  shift - positive parameter which depends on the time integration method
-.  V0 - time derivative of the state vector at t0
-.  t1 - time at step/stage being solved
-.  U1 - state vector at t1
-.  t0 - time at current step
-.  U0 - state vector at t0
-.  J - Jacobian matrix
--  ctx - [optional] user-defined context for evaluation routine
-
-   Level: normal
-
-.keywords: IGA, options
-@*/
-PetscErrorCode IGASetUserIEJacobian(IGA iga,IGAUserIEJacobian IEJacobian,void *IEJacCtx)
-{
-  PetscFunctionBegin;
-  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
-  if (IEJacobian) iga->userops->IEJacobian = IEJacobian;
-  if (IEJacCtx)   iga->userops->IEJacCtx   = IEJacCtx;
-  PetscFunctionReturn(0);
-}

File src/petigaaxis.c

   PetscFunctionBegin;
   PetscValidPointer(base,1);
   PetscValidPointer(axis,2);
-  if(base == axis) PetscFunctionReturn(0);
+  if (base == axis) PetscFunctionReturn(0);
 
   axis->periodic = base->periodic;
 

File src/petigabound.c

-#include "petiga.h"
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGABoundaryCreate"
-PetscErrorCode IGABoundaryCreate(IGABoundary *boundary)
-{
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidPointer(boundary,1);
-  ierr = PetscNew(struct _n_IGABoundary,boundary);CHKERRQ(ierr);
-  (*boundary)->refct = 1;
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__
-#define __FUNCT__ "IGABoundaryDestroy"
-PetscErrorCode IGABoundaryDestroy(IGABoundary *_boundary)
-{
-  IGABoundary    boundary;
-  PetscErrorCode ierr;
-  PetscFunctionBegin;
-  PetscValidPointer(_boundary,1);
-  boundary = *_boundary; *_boundary = 0;
-  if (!boundary) PetscFunctionReturn(0);
-  if (--boundary->refct > 0) PetscFunctionReturn(0);
-  ierr = IGABoundaryReset(boundary);CHKERRQ(ierr);
-  ierr = PetscFree(boundary);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef  __FUNCT__