Commits

Lisandro Dalcin  committed ef18631

Add support for TS problems with both implicit and explicit terms

  • Participants
  • Parent commits 930a7c5

Comments (0)

Files changed (6)

File demo/PatternFormation.c

 
 */
 
+#define EXPLICIT 1
 typedef struct {
+  PetscBool IMPLICIT;
   PetscReal delta;
   PetscReal D1,D2;
   PetscReal alpha;
 #undef  __FUNCT__
 #define __FUNCT__ "Function"
 PetscErrorCode Function(IGAPoint point,PetscReal dt,PetscReal shift,
-                        PetscReal t,const PetscScalar *V,const PetscScalar *U,
+                        PetscReal t,const PetscScalar *V,const PetscScalar *U1,const PetscScalar *U0,
                         PetscScalar *Re,void *ctx)
 {
   AppCtx *user = (AppCtx *)ctx;
-
-  PetscScalar uv_t[2],uv_0[2],uv_1[2][2];
-  IGAPointInterpolate(point,0,V,&uv_t[0]);
-  IGAPointInterpolate(point,0,U,&uv_0[0]);
-  IGAPointInterpolate(point,1,U,&uv_1[0][0]);
-  PetscScalar u_t = uv_t[0],    v_t = uv_t[1];
-  PetscScalar u   = uv_0[0],    v   = uv_0[1];
-  PetscScalar u_x = uv_1[0][0], v_x = uv_1[1][0];
-  PetscScalar u_y = uv_1[0][1], v_y = uv_1[1][1];
+  PetscBool IMPLICIT = user->IMPLICIT;
 
   PetscReal delta = user->delta;
   PetscReal D1    = user->D1;
   PetscReal tau1  = user->tau1;
   PetscReal tau2  = user->tau2;
 
+  PetscScalar uv_t[2],uv_0[2],uv_1[2][2];
+  IGAPointInterpolate(point,0,V,&uv_t[0]);
+  if (IMPLICIT)
+    IGAPointInterpolate(point,0,U1,&uv_0[0]);
+  else
+    IGAPointInterpolate(point,0,U0,&uv_0[0]);
+  IGAPointInterpolate(point,1,U1,&uv_1[0][0]);
+  PetscScalar u_t = uv_t[0],    v_t = uv_t[1];
+  PetscScalar u   = uv_0[0],    v   = uv_0[1];
+  PetscScalar u_x = uv_1[0][0], v_x = uv_1[1][0];
+  PetscScalar u_y = uv_1[0][1], v_y = uv_1[1][1];
+
   PetscReal f = alpha*u*(1-tau1*v*v) + v*(1-tau2*u);
   PetscReal g = beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);
   PetscReal (*N0)    = point->shape[0];
 #undef  __FUNCT__
 #define __FUNCT__ "Jacobian"
 PetscErrorCode Jacobian(IGAPoint point,PetscReal dt,PetscReal shift,
-                        PetscReal t,const PetscScalar *V,const PetscScalar *U,
+                        PetscReal t,const PetscScalar *V,const PetscScalar *U1, const PetscScalar *U0,
                         PetscScalar *Ke,void *ctx)
 {
   AppCtx *user = (AppCtx *)ctx;
-
-  PetscScalar uv_0[2];
-  IGAPointInterpolate(point,0,U,&uv_0[0]);
-  PetscScalar u = uv_0[0], v = uv_0[1];
+  PetscBool IMPLICIT = user->IMPLICIT;
 
   PetscReal delta = user->delta;
   PetscReal D1    = user->D1;
   PetscReal tau1  = user->tau1;
   PetscReal tau2  = user->tau2;
 
-  PetscReal f_u = alpha*(1-tau1*v*v) - tau2*v;
-  PetscReal f_v = -2*alpha*tau1*u*v + (1-tau2*u);
-  PetscReal g_u = alpha*tau1*v*v + (gamma+tau2*v);
-  PetscReal g_v = (beta+2*alpha*tau1*u*v) + tau2*u;
+  PetscReal f_u,f_v,g_u,g_v;
+  if (IMPLICIT) {
+    PetscScalar uv_0[2];
+    IGAPointInterpolate(point,0,U1,&uv_0[0]);
+    PetscScalar u = uv_0[0], v = uv_0[1];
+    f_u = alpha*(1-tau1*v*v) - tau2*v;
+    f_v = -2*alpha*tau1*u*v + (1-tau2*u);
+    g_u = alpha*tau1*v*v + (gamma+tau2*v);
+    g_v = (beta+2*alpha*tau1*u*v) + tau2*u;
+  }
 
   PetscReal (*N0)    = point->shape[0];
   PetscReal (*N1)[2] = (PetscReal (*)[2]) point->shape[1];
       PetscScalar Kab[2][2] = {{0,0},{0,0}};
       Kab[0][0] = shift*Na*Nb + delta*D1*(Na_x*Nb_x + Na_y*Nb_y);
       Kab[1][1] = shift*Na*Nb + delta*D2*(Na_x*Nb_x + Na_y*Nb_y);
-      Kab[0][0] -= Na*f_u*Nb; Kab[0][1] -= Na*f_v*Nb;
-      Kab[1][0] -= Na*g_u*Nb; Kab[1][1] -= Na*g_v*Nb;
-      for (i=0;i<2;i++)
-        for (j=0;j<2;j++)
-          K[a][i][b][j] += Kab[i][j];
+      if (IMPLICIT) {
+        Kab[0][0] -= Na*f_u*Nb; Kab[0][1] -= Na*f_v*Nb;
+        Kab[1][0] -= Na*g_u*Nb; Kab[1][1] -= Na*g_v*Nb;
+        for (i=0;i<2;i++)
+          for (j=0;j<2;j++)
+            K[a][i][b][j] += Kab[i][j];
+      } else {
+        K[a][0][b][0] += Kab[0][0];
+        K[a][1][b][1] += Kab[1][1];
+      }
     }
   }
   return 0;
   user.tau1  =  0.020;
   user.tau2  =  0.200;
 
+  user.IMPLICIT = PETSC_FALSE;
+
   PetscInt i;
   PetscInt dim = 2;
   PetscInt dof = 2;
   PetscInt p[2] = { 2, 2}, np = 2;
   PetscInt C[2] = {-1,-1}, nC = 2;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","PatternFormation Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-implicit","Treat all terms implicitly",__FILE__,user.IMPLICIT,&user.IMPLICIT,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsIntArray("-N","number of elements",     __FILE__,N,&nN,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsIntArray("-p","polynomial order",       __FILE__,p,&np,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsIntArray("-C","global continuity order",__FILE__,C,&nC,PETSC_NULL);CHKERRQ(ierr);
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  ierr = IGASetUserIFunction(iga,Function,&user);CHKERRQ(ierr);
-  ierr = IGASetUserIJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
+  ierr = IGASetUserIEFunction(iga,Function,&user);CHKERRQ(ierr);
+  ierr = IGASetUserIEJacobian(iga,Jacobian,&user);CHKERRQ(ierr);
 
   PetscReal h  = PetscMin(2.0/N[0],2.0/N[1]);
   PetscReal dt = h/user.delta/15.0;
   ierr = TSSetDuration(ts,120,10000.0);CHKERRQ(ierr);
   ierr = TSSetTime(ts,0.0);CHKERRQ(ierr);
   ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
+  if (!user.IMPLICIT) {
+    ierr = TSSetProblemType(ts,TS_LINEAR);CHKERRQ(ierr);
+  }
   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
 
   PetscReal t; Vec U;

File demo/makefile

 TARGETS = InputOutput \
-	  L2Proj1D L2Proj2D \
+          L2Proj1D L2Proj2D \
           Poisson1D Poisson2D Poisson3D \
           Laplace \
           Poisson \
 	-@${MPIEXEC} -n 1 ./CahnHilliard2D ${OPTS} -ts_max_steps 2
 runex4_4:
 	-@${MPIEXEC} -n 4 ./CahnHilliard2D ${OPTS} -ts_max_steps 2
+runex5a_1:
+	-@${MPIEXEC} -n 1 ./PatternFormation ${OPTS} -ts_max_steps 2
+runex5a_4:
+	-@${MPIEXEC} -n 4 ./PatternFormation ${OPTS} -ts_max_steps 2
+runex5b_1:
+	-@${MPIEXEC} -n 1 ./PatternFormation ${OPTS} -ts_max_steps 2 -implicit
+runex5b_4:
+	-@${MPIEXEC} -n 4 ./PatternFormation ${OPTS} -ts_max_steps 2 -implicit
 
 InputOutput := InputOutput.PETSc runex0a_1 runex0a_2 runex0a_3 runex0a_4 runex0a_8 runex0a.rm InputOutput.rm
 L2Proj1D := L2Proj1D.PETSc runex1a_1 runex1a_4 L2Proj1D.rm
 Poisson2D := Poisson2D.PETSc runex2b_1 runex2b_4 Poisson2D.rm
 Poisson3D := Poisson3D.PETSc runex2c_1 runex2c_4 Poisson3D.rm
 CahnHilliard2D := CahnHilliard2D.PETSc runex4_1 runex4_4 CahnHilliard2D.rm
+PatternFormation := PatternFormation.PETSc runex5a_1 runex5a_4 runex5b_1 runex5b_4 PatternFormation.rm
 
 L2Proj       := $(L2Proj1D) $(L2Proj2D) $(L2Proj3D)
 Poisson      := $(Poisson1D) $(Poisson2D) $(Poisson3D)
 CahnHilliard := $(CahnHilliard2D) $(CahnHilliard3D)
 
-TESTEXAMPLES_C := $(InputOutput) $(L2Proj) $(Poisson) $(CahnHilliard)
+TESTEXAMPLES_C := $(InputOutput) $(L2Proj) $(Poisson) $(CahnHilliard) $(PatternFormation)
 TESTEXAMPLES_F :=
 TESTEXAMPLES_FORTRAN:=$(TESTEXAMPLES_F)
 testexamples:
 	-@${OMAKE} tree ACTION=testexamples_C PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR}
 testfortran:
 	-@if [ "${FC}" != "" ]; then \
-            ${OMAKE} tree ACTION=testexamples_Fortran PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR}; \
-          fi
+	    ${OMAKE} tree ACTION=testexamples_Fortran PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR}; \
+	  fi
 
 include ${PETIGA_DIR}/conf/petigatest

File include/petiga.h

 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,PetscReal a,
-                                            PetscReal t,const PetscScalar *V,const PetscScalar *U,PetscScalar *F,
-                                            void *ctx);
+                                            PetscReal t,const PetscScalar *V,const PetscScalar *U,
+                                            PetscScalar *F,void *ctx);
 typedef PetscErrorCode (*IGAUserIJacobian) (IGAPoint point,PetscReal dt,PetscReal a,
-                                            PetscReal t,const PetscScalar *V,const PetscScalar *U,PetscScalar *J,
-                                            void *ctx);
+                                            PetscReal t,const PetscScalar *V,const PetscScalar *U,
+                                            PetscScalar *J,void *ctx);
+typedef PetscErrorCode (*IGAUserIEFunction)(IGAPoint point,PetscReal dt,PetscReal a,
+                                            PetscReal t,const PetscScalar *V,const PetscScalar *U,const PetscScalar *U0,
+                                            PetscScalar *F,void *ctx);
+typedef PetscErrorCode (*IGAUserIEJacobian)(IGAPoint point,PetscReal dt,PetscReal a,
+                                            PetscReal t,const PetscScalar *V,const PetscScalar *U,const PetscScalar *U0,
+                                            PetscScalar *J,void *ctx);
 
 typedef struct _IGAUserOps *IGAUserOps;
 struct _IGAUserOps {
   void              *IFunCtx;
   IGAUserIJacobian  IJacobian;
   void              *IJacCtx;
+  /**/
+  IGAUserIEFunction IEFunction;
+  void              *IEFunCtx;
+  IGAUserIEJacobian IEJacobian;
+  void              *IEJacCtx;
 };
 
 /* ---------------------------------------------------------------- */
 
 extern PetscErrorCode IGAGetElement(IGA iga,IGAElement *element);
 
-extern PetscErrorCode IGASetUserSystem   (IGA iga,IGAUserSystem    System,   void *SysCtx);
-extern PetscErrorCode IGASetUserFunction (IGA iga,IGAUserFunction  Function, void *FunCtx);
-extern PetscErrorCode IGASetUserJacobian (IGA iga,IGAUserJacobian  Jacobian, void *JacCtx);
-extern PetscErrorCode IGASetUserIFunction(IGA iga,IGAUserIFunction IFunction,void *FunCtx);
-extern PetscErrorCode IGASetUserIJacobian(IGA iga,IGAUserIJacobian IJacobian,void *JacCtx);
+extern PetscErrorCode IGASetUserSystem    (IGA iga,IGAUserSystem     System,    void *SysCtx);
+extern PetscErrorCode IGASetUserFunction  (IGA iga,IGAUserFunction   Function,  void *FunCtx);
+extern PetscErrorCode IGASetUserJacobian  (IGA iga,IGAUserJacobian   Jacobian,  void *JacCtx);
+extern PetscErrorCode IGASetUserIFunction (IGA iga,IGAUserIFunction  IFunction, void *FunCtx);
+extern PetscErrorCode IGASetUserIJacobian (IGA iga,IGAUserIJacobian  IJacobian, void *JacCtx);
+extern PetscErrorCode IGASetUserIEFunction(IGA iga,IGAUserIEFunction IEFunction,void *FunCtx);
+extern PetscErrorCode IGASetUserIEJacobian(IGA iga,IGAUserIEJacobian IEJacobian,void *JacCtx);
 
 /* ---------------------------------------------------------------- */
 
 extern PetscErrorCode IGAFormIJacobian(IGA iga,PetscReal dt,PetscReal a,
                                        PetscReal t,Vec V,Vec U,Mat J,
                                        IGAUserIJacobian IJacobian,void *ctx);
+extern PetscErrorCode IGAFormIEFunction(IGA iga,PetscReal dt,PetscReal a,
+                                        PetscReal t,Vec V,Vec U,Vec U0,Vec F,
+                                        IGAUserIEFunction IEFunction,void *ctx);
+extern PetscErrorCode IGAFormIEJacobian(IGA iga,PetscReal dt,PetscReal a,
+                                        PetscReal t,Vec V,Vec U,Vec U0,Mat J,
+                                        IGAUserIEJacobian IEJacobian,void *ctx);
 
 /* ---------------------------------------------------------------- */
 

File src/petiga.c

   PetscFunctionReturn(0);
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "IGASetUserIEFunction"
+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"
+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);
+}
 
 #if PETSC_VERSION_(3,2,0)
 #include "private/dmimpl.h"

File src/petigaelem.c

   PetscFunctionReturn(0);
 }
 
-
 #undef  __FUNCT__
 #define __FUNCT__ "IGAElementAssembleVec"
 PetscErrorCode IGAElementAssembleVec(IGAElement element,PetscScalar F[],Vec vec)

File src/petigats.c

 #undef  __FUNCT__
 #define __FUNCT__ "IGAFormIFunction"
 PetscErrorCode IGAFormIFunction(IGA iga,PetscReal dt,PetscReal shift,
-                                PetscReal t,Vec vecV,Vec vecU,Vec vecF,
+                                PetscReal t,Vec vecV,Vec vecU,
+                                Vec vecF,
                                 IGAUserIFunction IFunction, void *ctx)
 {
   Vec               localV;
 #undef  __FUNCT__
 #define __FUNCT__ "IGAFormIJacobian"
 PetscErrorCode IGAFormIJacobian(IGA iga,PetscReal dt,PetscReal shift,
-                                PetscReal t,Vec vecV,Vec vecU,Mat matJ,
+                                PetscReal t,Vec vecV,Vec vecU,
+                                Mat matJ,
                                 IGAUserIJacobian IJacobian,void *ctx)
 {
   Vec               localV;
   ierr = IGAElementBegin(element);CHKERRQ(ierr);
   while (IGAElementNext(element)) {
     PetscScalar *V, *U, *J;
+    ierr = IGAElementGetWorkVec(element,&V);CHKERRQ(ierr);
     ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
-    ierr = IGAElementGetWorkVec(element,&V);CHKERRQ(ierr);
     ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
     /* */
     ierr = IGAElementGetValues(element,arrayV,V);CHKERRQ(ierr);
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAFormIEFunction"
+PetscErrorCode IGAFormIEFunction(IGA iga,PetscReal dt,PetscReal shift,
+                                 PetscReal t,Vec vecV,Vec vecU,Vec vecU0,
+                                 Vec vecF,
+                                 IGAUserIEFunction IEFunction, void *ctx)
+{
+  Vec               localV;
+  Vec               localU;
+  Vec               localU0;
+  const PetscScalar *arrayV;
+  const PetscScalar *arrayU;
+  const PetscScalar *arrayU0;
+  IGAElement        element;
+  IGAPoint          point;
+  PetscErrorCode    ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(vecV,VEC_CLASSID,5);
+  PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
+  PetscValidHeaderSpecific(vecU0,VEC_CLASSID,7);
+  PetscValidHeaderSpecific(vecF,VEC_CLASSID,8);
+  IGACheckSetUp(iga,1);
+
+  /* Clear global vector F */
+  ierr = VecZeroEntries(vecF);CHKERRQ(ierr);
+
+  /* Get local vector and array V */
+  ierr = IGAGetLocalVec(iga,&localV);CHKERRQ(ierr);
+  ierr = IGAGlobalToLocal(iga,vecV,localV);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(localV,&arrayV);CHKERRQ(ierr);
+  /* Get local vector and array U */
+  ierr = IGAGetLocalVec(iga,&localU);CHKERRQ(ierr);
+  ierr = IGAGlobalToLocal(iga,vecU,localU);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(localU,&arrayU);CHKERRQ(ierr);
+  /* Get local vector and array U0 */
+  ierr = IGAGetLocalVec(iga,&localU0);CHKERRQ(ierr);
+  ierr = IGAGlobalToLocal(iga,vecU0,localU0);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(localU0,&arrayU0);CHKERRQ(ierr);
+
+  /* Element loop */
+  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
+  ierr = IGAElementBegin(element);CHKERRQ(ierr);
+  while (IGAElementNext(element)) {
+    PetscScalar *V, *U, *U0, *F;
+    ierr = IGAElementGetWorkVec(element,&V);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&U0);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&F);CHKERRQ(ierr);
+    /* */
+    ierr = IGAElementGetValues(element,arrayV,V);CHKERRQ(ierr);
+    ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
+    ierr = IGAElementGetValues(element,arrayU0,U0);CHKERRQ(ierr);
+    ierr = IGAElementFixValues(element,U0);CHKERRQ(ierr);
+    ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
+    /* Quadrature loop */
+    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
+    ierr = IGAPointBegin(point);CHKERRQ(ierr);
+    while (IGAPointNext(point)) {
+      PetscScalar *R;
+      ierr = IGAPointGetWorkVec(point,&R);CHKERRQ(ierr);
+      ierr = IEFunction(point,dt,shift,t,V,U,U0,R,ctx);CHKERRQ(ierr);
+      ierr = IGAPointAddVec(point,R,F);CHKERRQ(ierr);
+    }
+    /* */
+    ierr = IGAElementFixFunction(element,F);CHKERRQ(ierr);
+    ierr = IGAElementAssembleVec(element,F,vecF);CHKERRQ(ierr);
+  }
+  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+
+  /* Restore array and local vector V */
+  ierr = VecRestoreArrayRead(localV,&arrayV);CHKERRQ(ierr);
+  ierr = IGARestoreLocalVec(iga,&localV);CHKERRQ(ierr);
+  /* Restore array and local vector U */
+  ierr = VecRestoreArrayRead(localU,&arrayU);CHKERRQ(ierr);
+  ierr = IGARestoreLocalVec(iga,&localU);CHKERRQ(ierr);
+  /* Restore array and local vector U0 */
+  ierr = VecRestoreArrayRead(localU0,&arrayU0);CHKERRQ(ierr);
+  ierr = IGARestoreLocalVec(iga,&localU0);CHKERRQ(ierr);
+
+  /* Assemble global vector F */
+  ierr = VecAssemblyBegin(vecF);CHKERRQ(ierr);
+  ierr = VecAssemblyEnd(vecF);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAFormIEJacobian"
+PetscErrorCode IGAFormIEJacobian(IGA iga,PetscReal dt,PetscReal shift,
+                                 PetscReal t,Vec vecV,Vec vecU,Vec vecU0,
+                                 Mat matJ,
+                                 IGAUserIEJacobian IEJacobian,void *ctx)
+{
+  Vec               localV;
+  Vec               localU;
+  Vec               localU0;
+  const PetscScalar *arrayV;
+  const PetscScalar *arrayU;
+  const PetscScalar *arrayU0;
+  IGAElement        element;
+  IGAPoint          point;
+  PetscErrorCode    ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
+  PetscValidHeaderSpecific(vecV,VEC_CLASSID,5);
+  PetscValidHeaderSpecific(vecU,VEC_CLASSID,6);
+  PetscValidHeaderSpecific(vecU0,VEC_CLASSID,7);
+  PetscValidHeaderSpecific(matJ,MAT_CLASSID,0);
+  IGACheckSetUp(iga,1);
+
+  /* Clear global matrix J*/
+  ierr = MatZeroEntries(matJ);CHKERRQ(ierr);
+
+  /* Get local vector and array V */
+  ierr = IGAGetLocalVec(iga,&localV);CHKERRQ(ierr);
+  ierr = IGAGlobalToLocal(iga,vecV,localV);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(localV,&arrayV);CHKERRQ(ierr);
+  /* Get local vector and array U */
+  ierr = IGAGetLocalVec(iga,&localU);CHKERRQ(ierr);
+  ierr = IGAGlobalToLocal(iga,vecU,localU);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(localU,&arrayU);CHKERRQ(ierr);
+  /* Get local vector and array U0 */
+  ierr = IGAGetLocalVec(iga,&localU0);CHKERRQ(ierr);
+  ierr = IGAGlobalToLocal(iga,vecU0,localU0);CHKERRQ(ierr);
+  ierr = VecGetArrayRead(localU0,&arrayU0);CHKERRQ(ierr);
+
+  /* Element Loop */
+  ierr = IGAGetElement(iga,&element);CHKERRQ(ierr);
+  ierr = IGAElementBegin(element);CHKERRQ(ierr);
+  while (IGAElementNext(element)) {
+    PetscScalar *V, *U, *U0, *J;
+    ierr = IGAElementGetWorkVec(element,&V);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&U);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkVec(element,&U0);CHKERRQ(ierr);
+    ierr = IGAElementGetWorkMat(element,&J);CHKERRQ(ierr);
+    /* */
+    ierr = IGAElementGetValues(element,arrayV,V);CHKERRQ(ierr);
+    ierr = IGAElementGetValues(element,arrayU,U);CHKERRQ(ierr);
+    ierr = IGAElementGetValues(element,arrayU0,U0);CHKERRQ(ierr);
+    ierr = IGAElementFixValues(element,U0);CHKERRQ(ierr);
+    ierr = IGAElementFixValues(element,U);CHKERRQ(ierr);
+    /* Quadrature loop */
+    ierr = IGAElementGetPoint(element,&point);CHKERRQ(ierr);
+    ierr = IGAPointBegin(point);CHKERRQ(ierr);
+    while (IGAPointNext(point)) {
+      PetscScalar *K;
+      ierr = IGAPointGetWorkMat(point,&K);CHKERRQ(ierr);
+      ierr = IEJacobian(point,dt,shift,t,V,U,U0,K,ctx);CHKERRQ(ierr);
+      ierr = IGAPointAddMat(point,K,J);CHKERRQ(ierr);
+    }
+    /* */
+    ierr = IGAElementFixJacobian(element,J);CHKERRQ(ierr);
+    ierr = IGAElementAssembleMat(element,J,matJ);CHKERRQ(ierr);
+  }
+  ierr = IGAElementEnd(element);CHKERRQ(ierr);
+
+  /* Restore array and local vector V */
+  ierr = VecRestoreArrayRead(localV,&arrayV);CHKERRQ(ierr);
+  ierr = IGARestoreLocalVec(iga,&localV);CHKERRQ(ierr);
+  /* Restore array and local vector U */
+  ierr = VecRestoreArrayRead(localU,&arrayU);CHKERRQ(ierr);
+  ierr = IGARestoreLocalVec(iga,&localU);CHKERRQ(ierr);
+  /* Restore array and local vector U0 */
+  ierr = VecRestoreArrayRead(localU0,&arrayU0);CHKERRQ(ierr);
+  ierr = IGARestoreLocalVec(iga,&localU0);CHKERRQ(ierr);
+
+  /* Assemble global matrix J*/
+  ierr = MatAssemblyBegin(matJ,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd  (matJ,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGATSFormIFunction"
 PetscErrorCode IGATSFormIFunction(TS ts,PetscReal t,Vec U,Vec V,Vec F,void *ctx)
 {
   PetscValidHeaderSpecific(V,VEC_CLASSID,4);
   PetscValidHeaderSpecific(F,VEC_CLASSID,5);
   PetscValidHeaderSpecific(iga,IGA_CLASSID,6);
-  if (!iga->userops->IFunction)
+  if (!iga->userops->IFunction && !iga->userops->IEFunction)
     SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call IGASetUserIFunction()");
   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
-  ierr = IGAFormIFunction(iga,dt,shift,t,V,U,F,
-                          iga->userops->IFunction,
-                          iga->userops->IFunCtx);CHKERRQ(ierr);
+  if (iga->userops->IEFunction) {
+    Vec U0;
+    ierr = TSGetSolution(ts,&U0);CHKERRQ(ierr);
+    ierr = IGAFormIEFunction(iga,dt,shift,t,V,U,U0,F,
+                             iga->userops->IEFunction,
+                             iga->userops->IEFunCtx);CHKERRQ(ierr);
+  } else {
+    ierr = IGAFormIFunction(iga,dt,shift,t,V,U,F,
+                            iga->userops->IFunction,
+                            iga->userops->IFunCtx);CHKERRQ(ierr);
+  }
   PetscFunctionReturn(0);
 }
 
   PetscValidHeaderSpecific(*P,MAT_CLASSID,7);
   PetscValidPointer(m,8);
   PetscValidHeaderSpecific(iga,IGA_CLASSID,9);
-  if (!iga->userops->IJacobian)
+  if (!iga->userops->IJacobian && !iga->userops->IEJacobian)
     SETERRQ(((PetscObject)ts)->comm,PETSC_ERR_USER,"Must call IGASetUserIJacobian()");
   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
-  ierr = IGAFormIJacobian(iga,dt,shift,t,V,U,*P,
-                          iga->userops->IJacobian,
-                          iga->userops->IJacCtx);CHKERRQ(ierr);
+  if (iga->userops->IEJacobian) {
+    Vec U0;
+    ierr = TSGetSolution(ts,&U0);CHKERRQ(ierr);
+    ierr = IGAFormIEJacobian(iga,dt,shift,t,V,U,U0,*P,
+                             iga->userops->IEJacobian,
+                             iga->userops->IEJacCtx);CHKERRQ(ierr);
+  } else {
+    ierr = IGAFormIJacobian(iga,dt,shift,t,V,U,*P,
+                            iga->userops->IJacobian,
+                            iga->userops->IJacCtx);CHKERRQ(ierr);
+  }
   if (*J != * P) {
     ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
     ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);