Commits

Nathan Collier committed fd5d1b6 Merge

branch merge

  • Participants
  • Parent commits e2ae03c, 5de1176

Comments (0)

Files changed (81)

   set (CMAKE_MODULE_PATH ${PetIGA_SOURCE_DIR}/conf/cmake)
   find_package (PETSc)
 else ()
-  find_path (PETSC_DIR include/petsc.h HINTS ENV PETSC_DIR DOC "PETSc Directory")
-  set (PETSC_ARCH $ENV{PETSC_ARCH})
+  find_path (PETSC_DIR include/petsc.h HINTS ENV PETSC_DIR PATHS $ENV{HOME}/petsc DOC "PETSc top-level directory")
+  set (PETSC_ARCH $ENV{PETSC_ARCH} CACHE STRING "PETSc configuration")
   find_path (PETSC_INCLUDE_DIR  petsc.h HINTS "${PETSC_DIR}" PATH_SUFFIXES include NO_DEFAULT_PATH)
-  find_path (PETSC_INCLUDE_CONF petscconf.h HINTS "${PETSC_DIR}" PATH_SUFFIXES "${PETSC_ARCH}/include" NO_DEFAULT_PATH)
+  find_path (PETSC_INCLUDE_CONF petscconf.h HINTS "${PETSC_DIR}" PATH_SUFFIXES "${PETSC_ARCH}/include" "include" NO_DEFAULT_PATH)
   mark_as_advanced (PETSC_INCLUDE_DIR PETSC_INCLUDE_CONF)
   set (PETSC_INCLUDES ${PETSC_INCLUDE_CONF} ${PETSC_INCLUDE_DIR} CACHE PATH "PETSc include paths" FORCE)
   set (PETSC_DEFINITIONS "-D__INSDIR__=" CACHE STRING "PETSc preprocesor definitions" FORCE)
-  find_path (PETSC_LIB_DIR NAMES "" HINTS "${PETSC_DIR}" PATH_SUFFIXES "${PETSC_ARCH}/lib" NO_DEFAULT_PATH)
-  find_library (PETSC_LIBRARY NAMES petsc HINTS "${PETSC_LIB_DIR}" NO_DEFAULT_PATH)
-  set (PETSC_LIBRARIES ${PETSC_LIBRARY} CACHE FILEPATH "PETSc library" FORCE)
+  mark_as_advanced (PETSC_DEFINITIONS)
+  find_library (PETSC_LIBRARIES NAMES petsc HINTS "${PETSC_DIR}" PATH_SUFFIXES "${PETSC_ARCH}/lib" "lib" NO_DEFAULT_PATH)
   include (${PETSC_DIR}/${PETSC_ARCH}/conf/PETScConfig.cmake)
-endif()
+  mark_as_advanced (PETSC_CLANGUAGE_Cxx)
+endif ()
 
 if (PETSC_CLANGUAGE_Cxx)
   enable_language (CXX)
                     "${PetIGA_SOURCE_DIR}/include")
 add_definitions (${PETSC_DEFINITIONS})
 
-set (CMAKE_ARCHIVE_OUTPUT_DIRECTORY "${PetIGA_BINARY_DIR}/lib" CACHE PATH "Output directory for PetIGA archives")
-set (CMAKE_LIBRARY_OUTPUT_DIRECTORY "${PetIGA_BINARY_DIR}/lib" CACHE PATH "Output directory for PetIGA libraries")
-set (CMAKE_Fortran_MODULE_DIRECTORY "${PetIGA_BINARY_DIR}/include" CACHE PATH "Output directory for fortran *.mod files")
+set (CMAKE_ARCHIVE_OUTPUT_DIRECTORY "${PetIGA_BINARY_DIR}/lib"     CACHE PATH "Output directory for PetIGA archives")
+set (CMAKE_LIBRARY_OUTPUT_DIRECTORY "${PetIGA_BINARY_DIR}/lib"     CACHE PATH "Output directory for PetIGA libraries")
+set (CMAKE_Fortran_MODULE_DIRECTORY "${PetIGA_BINARY_DIR}/include" CACHE PATH "Output directory for Fortran modules")
 mark_as_advanced (CMAKE_ARCHIVE_OUTPUT_DIRECTORY CMAKE_LIBRARY_OUTPUT_DIRECTORY CMAKE_Fortran_MODULE_DIRECTORY)
 set (CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
 set (CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
 file (GLOB PetIGA_SOURCES_C RELATIVE ${PetIGA_SOURCE_DIR} ${PetIGA_SOURCE_DIR}/src/*.c)
 file (GLOB PetIGA_SOURCES_F RELATIVE ${PetIGA_SOURCE_DIR} ${PetIGA_SOURCE_DIR}/src/*.[Ff]90)
 set  (PetIGA_SOURCES_ALL ${PetIGA_SOURCES_C} ${PetIGA_SOURCES_F})
-
-add_library (petiga ${PetIGA_SOURCES_ALL})
-target_link_libraries (petiga ${PETSC_LIBRARIES} ${PETSC_PACKAGE_LIBS})
-
 if (PETSC_CLANGUAGE_Cxx)
   foreach (file ${PetIGA_SOURCES_C})
     set_source_files_properties(${file} PROPERTIES LANGUAGE CXX)
   endforeach ()
 endif ()
 
+add_library (petiga ${PetIGA_SOURCES_ALL})
+target_link_libraries (petiga ${PETSC_LIBRARIES} ${PETSC_PACKAGE_LIBS})
+
 #set (BUILD_SHARED_LIBS NO)
 
 #if (BUILD_SHARED_LIBS)
 set appropriate values for ``PETSC_DIR`` and ``PETSC_ARCH`` in your
 environment::
 
-  $ export PETSC_DIR=/home/user/petsc-3.3-p3
+  $ export PETSC_DIR=/home/user/petsc-3.4.0
   $ export PETSC_ARCH=arch-linux2-c-debug
 
 Clone the `Mercurial <http://mercurial.selenic.com/>`_ repository

demo/AdaptiveL2Projection.c

+/*
+  This code computes the L2 projection of a function defined in
+  Function to a B-spline space. The space is adaptively refined by a
+  brute force strategy: which knot do we insert such that the global
+  L2 error is reduced the most?
+
+  keywords: scalar, linear
+ */
 #include "petiga.h"
 
 #define SQ(A) (A)*(A)
 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen;
-  IGAPointGetSizes(p,0,&nen,0);
+  PetscInt nen = p->nen;
   
   PetscReal x[3] = {0,0,0};
   IGAPointFormPoint(p,x);
 #define __FUNCT__ "ComputeError"
 PetscErrorCode ComputeError(PetscInt dim,PetscInt p,PetscInt C,PetscReal (*U)[MAX_BREAKS],PetscInt *N,PetscReal *error)
 {
-  
   PetscErrorCode  ierr;
   PetscInt i;
   IGA iga;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
   PetscReal errors[2];
-  ierr  = IGAFormScalar(iga,x,2,&errors[0],Error,PETSC_NULL);CHKERRQ(ierr);
+  ierr  = IGAFormScalar(iga,x,2,&errors[0],Error,NULL);CHKERRQ(ierr);
   errors[0] = PetscSqrtReal(PetscRealPart(errors[0]));
   errors[1] = PetscSqrtReal(PetscRealPart(errors[1]));
   *error = errors[0]/errors[1];
   PetscInt i,j;
   PetscInt dim = 2, p = 2, C = 1;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-p","polynomial order",       __FILE__,p,&p,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-C","global continuity order",__FILE__,C,&C,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-p","polynomial order",       __FILE__,p,&p,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-C","global continuity order",__FILE__,C,&C,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   // initial breaks (2x2 mesh)
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  ierr = IGASetUserSystem(iga,System,NULL);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
 
   KSP ksp;

demo/AdvectionDiffusion.c

+/*
+  This code solves the advection-diffusion problem where the advection
+  is constant and skew to the mesh. It is independent of the spatial
+  dimension which must be specified in the commandline options via
+  -iga_dim. The strength of advection can be controlled through the
+  option -Pe, representing the Peclet number. For example,
+
+  ./AdvectionDiffusion -iga_dim 1 -draw -draw_pause -1 -Pe 10
+
+  keywords: steady, scalar, linear, educational, dimension independent
+ */
 #include "petiga.h"
 
 typedef struct {
   PetscReal adv[3];
 } AppCtx;
 
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
   AppCtx *user = (AppCtx *)ctx;
 
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
-
-  const PetscReal *N;
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,0,&N);
-  IGAPointGetShapeFuns(p,1,&N1);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+  const PetscReal *N0        = (typeof(N0)) p->shape[0];
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
 
   PetscInt a,b,i;
   for (a=0; a<nen; a++) {
     for (b=0; b<nen; b++) {
-      PetscScalar diffusion = 0.0;
+      PetscScalar diffusion = DOT(dim,N1[a],N1[b]);
       PetscScalar advection = 0.0;
-      for (i=0; i<dim; i++){
-        diffusion += N1[a*dim+i]*N1[b*dim+i];
-	advection += user->adv[i]*N1[b*dim+i];
-      }
-      advection *= N[a];
-      K[a*nen+b] = diffusion + advection; 
+      for (i=0; i<dim; i++)
+        advection += user->adv[i]*N1[b][i];
+      advection *= N0[a];
+      K[a*nen+b] = diffusion + advection;
     }
     F[a] = 0.0;
   }
   PetscErrorCode ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
-  AppCtx user;
-  PetscInt  i;
-  PetscInt  dim = 3;
-  PetscInt  dof = 1;
-  PetscInt  N[3] = {16,16,16};
-  PetscInt  p[3] = { 2, 2, 2};
-  PetscInt  C[3] = {-1,-1,-1};
-  PetscInt  n1=3, n2=3, n3=3;
-  PetscReal Pe=1.0;
+  AppCtx    user;
+  PetscInt  i,dim;
+  PetscReal Pe = 1.0;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","AdvectionDiffusion Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
-  n1 = n2 = n3 = dim;
-  ierr = PetscOptionsIntArray ("-N","number of elements",     __FILE__,N,&n1,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray ("-p","polynomial order",       __FILE__,p,&n2,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray ("-C","global continuity order",__FILE__,C,&n3,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsReal     ("-Pe","Peclet number",__FILE__,Pe,&Pe,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-Pe","Peclet number",__FILE__,Pe,&Pe,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
-  if (n1<3) N[2] = N[0]; if (n1<2) N[1] = N[0];
-  if (n2<3) p[2] = p[0]; if (n2<2) p[1] = p[0];
-  if (n3<3) C[2] = C[0]; if (n3<2) C[1] = C[0];
-  for (i=0; i<dim; i++)  if (C[i] ==-1) C[i] = p[i] - 1;
-  for (i=0; i<dim; i++) user.adv[i] = Pe*sqrt(dim)/dim;
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
-  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
-  ierr = IGASetDof(iga,dof);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  IGABoundary bnd;
   for (i=0; i<dim; i++) {
-    IGAAxis axis;
-    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);
-    IGABoundary bnd;
     ierr = IGAGetBoundary(iga,i,0,&bnd);CHKERRQ(ierr);
     ierr = IGABoundarySetValue(bnd,0,1.0);CHKERRQ(ierr);
     ierr = IGAGetBoundary(iga,i,1,&bnd);CHKERRQ(ierr);
     ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
   }
-  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  // advection is skew to the mesh
+  for (i=0; i<dim; i++) user.adv[i] = Pe*sqrt(dim)/dim;
 
   Mat A;
   Vec x,b;

demo/BoundaryIntegral.c

+/*
+  This code solves the Laplace problem where the boundary conditions
+  can be changed from Neumann to Dirichlet via the commandline. While
+  its primary use is in regression tests for PetIGA, it also
+  demonstrates how boundary integrals may be performed to enforce
+  things like Neumann conditions.
+
+  keywords: steady, scalar, linear, testing, dimension independent,
+  boundary integrals
+ */
 #include "petiga.h"
 
+typedef struct {
+  PetscInt dir;
+  PetscInt side;
+} AppCtx;
+
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
-
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
-
-  PetscInt a,b,i;
+  PetscInt dim = p->dim;
+  PetscInt nen = p->nen;
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
+  PetscInt a,b;
   for (a=0; a<nen; a++) {
-    for (b=0; b<nen; b++) {
-      PetscScalar Kab = 0.0;
-      for (i=0; i<dim; i++)
-        Kab += N1[a*dim+i]*N1[b*dim+i];
-      K[a*nen+b] = Kab;
-    }
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
+    F[a] = 0.0;
   }
   return 0;
 }
 
-
-
 #undef  __FUNCT__
 #define __FUNCT__ "Neumann"
 PetscErrorCode Neumann(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscReal *N0 = p->shape[0];
-  PetscInt a,nen=p->nen;
-  for (a=0; a<nen; a++) {
-    PetscReal Na   = N0[a];
-    F[a] = Na * 1.0;
-  }
+  PetscInt nen = p->nen;
+  const PetscReal *N0 = (typeof(N0)) p->shape[0];
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    F[a] = N0[a] * 1.0;
   return 0;
 }
 
-typedef struct { 
-  PetscInt dir;
-  PetscInt side;
-} AppCtx;
+PETSC_STATIC_INLINE
+PetscReal DEL2(PetscInt dim,const PetscReal a[dim][dim])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i][i];
+  return s;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "SystemCollocation"
+PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+  const PetscReal (*N2)[dim][dim] = (typeof(N2)) p->shape[2];
+
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = -DEL2(dim,N2[a]);
+  F[0] = 0.0;
+
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "DirichletCollocation"
+PetscErrorCode DirichletCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen = p->nen;
+  const PetscReal *N0 = (typeof(N0)) p->shape[0];
+
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = N0[a];
+  F[0] = 1.0;
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "NeumannCollocation"
+PetscErrorCode NeumannCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
+  const PetscReal *normal    = p->normal;
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = DOT(dim,N1[a],normal);
+  F[0] = 1.0;
+  return 0;
+}
+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__ "Error"
   IGAPointFormValue(p,U,&u);
   PetscReal x;
   if (user->side == 0)
-    x = 1 - p->point[user->dir];
+    x = 1 - p->point[user->dir] + 1;
   else
-    x = p->point[user->dir];
+    x = p->point[user->dir] + 1;
   PetscReal e = u - x;
   S[0] = e*e;
   return 0;
   user.dir  = 0;
   user.side = 1;
 
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-dir", "direction",__FILE__,user.dir, &user.dir, PETSC_NULL);CHKERRQ(ierr); 
-  ierr = PetscOptionsInt("-side","side",     __FILE__,user.side,&user.side,PETSC_NULL);CHKERRQ(ierr); 
+  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 ("-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);
+  ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  if (iga->dim < 1) {ierr = IGASetDim(iga,2);CHKERRQ(ierr);}
 
-  IGABoundary bnd;
-  PetscInt d = !user.side; 
+  PetscInt dim;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  PetscInt d =  !user.side;
   PetscInt n = !!user.side;
-  ierr = IGAGetBoundary(iga,user.dir,d,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-  ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr); 
-  ierr = IGABoundarySetUserSystem(bnd,Neumann,PETSC_NULL);CHKERRQ(ierr);
-
-  ierr = IGASetFromOptions(iga);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);
+  } 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);
+  }
   ierr = IGASetUp(iga);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,PETSC_NULL);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 = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
   ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
-  PetscInt dim;
-  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
-  if (dim <= 2) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
-
   PetscScalar error = 0;
   ierr = IGAFormScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
-  ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);
-  
+
+  if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);}
+  if (check_error) {if (error>1e-4) SETERRQ1(PETSC_COMM_WORLD,1,"L2 error=%G\n",error);}
+  if (draw&&dim<3) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);
   ierr = VecDestroy(&b);CHKERRQ(ierr);
   ierr = IGADestroy(&iga);CHKERRQ(ierr);
 
-  PetscBool flag = PETSC_FALSE;
-  PetscReal secs = -1;
-  ierr = PetscOptionsHasName(0,"-sleep",&flag);CHKERRQ(ierr);
-  ierr = PetscOptionsGetReal(0,"-sleep",&secs,0);CHKERRQ(ierr);
-  if (flag) {ierr = PetscSleep(secs);CHKERRQ(ierr);}
-
   ierr = PetscFinalize();CHKERRQ(ierr);
   return 0;
 }
+/*
+   This code solves the steady and unsteady Bratu equation. It also
+   demonstrates how the user-specified routines, here the Function and
+   Jacobian routines, can be implemented in Fortran (see BratuFJ.F90)
+   yet called from PetIGA in C.
+
+   keywords: steady, transient, scalar, implicit, nonlinear, testing,
+   dimension independent, collocation, fortran
+*/
 #include "petiga.h"
 
-typedef struct { 
+typedef struct {
   PetscReal lambda;
 } AppCtx;
 
 
   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
-  
-  
+
+
   PetscBool steady = PETSC_TRUE;
   PetscReal lambda = 6.80;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Bratu Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-lambda","Bratu parameter",__FILE__,lambda,&lambda,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-steady","Steady problem",__FILE__,steady,&steady,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-lambda","Bratu parameter",__FILE__,lambda,&lambda,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-steady","Steady problem",__FILE__,steady,&steady,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
-
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-  IGABoundary bnd;
   PetscInt dir,side;
   for (dir=0; dir<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,0,0.0);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,field,value);CHKERRQ(ierr);
     }
   }
-
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-
   PetscInt dim;
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   if (dim < 1) {ierr = IGASetDim(iga,dim=2);CHKERRQ(ierr);}
   if (steady) {
     SNES snes;
     ierr = IGACreateSNES(iga,&snes);CHKERRQ(ierr);
+    ierr = SNESSetTolerances(snes,PETSC_DEFAULT,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
     ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
     ierr = SNESSolve(snes,0,x);CHKERRQ(ierr);
     ierr = SNESDestroy(&snes);CHKERRQ(ierr);
   } else {
     TS ts;
+    SNES snes;
     ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
     ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr);
     ierr = TSSetDuration(ts,10000,0.1);CHKERRQ(ierr);
     ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr);
+    ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
+    ierr = SNESSetTolerances(snes,PETSC_DEFAULT,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
     ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
 #if PETSC_VERSION_LE(3,3,0)
-    ierr = TSSolve(ts,x,PETSC_NULL);CHKERRQ(ierr);
+    ierr = TSSolve(ts,x,NULL);CHKERRQ(ierr);
 #else
     ierr = TSSolve(ts,x);CHKERRQ(ierr);
 #endif

demo/CahnHilliard2D.c

   PetscBool monitor = PETSC_FALSE; 
   char initial[PETSC_MAX_PATH_LEN] = {0};
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","CahnHilliard2D Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-N","number of elements (along one dimension)",__FILE__,N,&N,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-p","polynomial order",__FILE__,p,&p,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-C","global continuity order",__FILE__,C,&C,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsString("-ch_initial","Load initial solution from file",__FILE__,initial,initial,sizeof(initial),PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-ch_output","Enable output files",__FILE__,output,&output,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-ch_monitor","Compute and show statistics of solution",__FILE__,monitor,&monitor,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-ch_cbar","Initial average concentration",__FILE__,user.cbar,&user.cbar,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-ch_alpha","Characteristic parameter",__FILE__,user.alpha,&user.alpha,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-N","number of elements (along one dimension)",__FILE__,N,&N,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-p","polynomial order",__FILE__,p,&p,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-C","global continuity order",__FILE__,C,&C,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsString("-ch_initial","Load initial solution from file",__FILE__,initial,initial,sizeof(initial),NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-ch_output","Enable output files",__FILE__,output,&output,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-ch_monitor","Compute and show statistics of solution",__FILE__,monitor,&monitor,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-ch_cbar","Initial average concentration",__FILE__,user.cbar,&user.cbar,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-ch_alpha","Characteristic parameter",__FILE__,user.alpha,&user.alpha,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
   if (C == PETSC_DECIDE) C = p-1;
 
 
   ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
   ierr = TSAlphaSetRadius(ts,0.5);CHKERRQ(ierr);
-  ierr = TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,PETSC_NULL);CHKERRQ(ierr); 
+  ierr = TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,NULL);CHKERRQ(ierr); 
 
   if (monitor) {
     user.iga = iga;
     PetscPrintf(PETSC_COMM_WORLD,"#Time        dt           Free Energy            Second moment          Third moment\n");
-    ierr = TSMonitorSet(ts,StatsMonitor,&user,PETSC_NULL);CHKERRQ(ierr);
+    ierr = TSMonitorSet(ts,StatsMonitor,&user,NULL);CHKERRQ(ierr);
   }
   if (output) {
-    ierr = TSMonitorSet(ts,OutputMonitor,&user,PETSC_NULL);CHKERRQ(ierr);
+    ierr = TSMonitorSet(ts,OutputMonitor,&user,NULL);CHKERRQ(ierr);
   }
   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
 
   ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
   ierr = FormInitialCondition(&user,iga,initial,U);CHKERRQ(ierr);
 #if PETSC_VERSION_LE(3,3,0)
-  ierr = TSSolve(ts,U,PETSC_NULL);CHKERRQ(ierr);
+  ierr = TSSolve(ts,U,NULL);CHKERRQ(ierr);
 #else
   ierr = TSSolve(ts,U);CHKERRQ(ierr);
 #endif

demo/CahnHilliard3D.c

   PetscBool monitor = PETSC_FALSE; 
   char initial[PETSC_MAX_PATH_LEN] = {0};
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","CahnHilliard3D Options","CH3D");CHKERRQ(ierr);
-  ierr = PetscOptionsString("-ch_initial","Load initial solution from file",__FILE__,initial,initial,sizeof(initial),PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-ch_output","Enable output files",__FILE__,output,&output,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-ch_monitor","Compute and show statistics of solution",__FILE__,monitor,&monitor,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-ch_cbar","Initial average concentration",__FILE__,user.cbar,&user.cbar,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-ch_alpha","Characteristic parameter",__FILE__,user.alpha,&user.alpha,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsString("-ch_initial","Load initial solution from file",__FILE__,initial,initial,sizeof(initial),NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-ch_output","Enable output files",__FILE__,output,&output,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-ch_monitor","Compute and show statistics of solution",__FILE__,monitor,&monitor,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-ch_cbar","Initial average concentration",__FILE__,user.cbar,&user.cbar,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-ch_alpha","Characteristic parameter",__FILE__,user.alpha,&user.alpha,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   IGA iga;
   ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
   ierr = TSAlphaSetRadius(ts,0.5);CHKERRQ(ierr);
 
-  user.E = PETSC_NULL;
+  user.E = NULL;
   ierr = TSAlphaSetAdapt(ts,CHAdapt,&user);CHKERRQ(ierr); 
 
   if (monitor) {
     user.iga = iga;
     PetscPrintf(PETSC_COMM_WORLD,"#Time        dt           Free Energy            Second moment          Third moment\n");
-    ierr = TSMonitorSet(ts,StatsMonitor,&user,PETSC_NULL);CHKERRQ(ierr);
+    ierr = TSMonitorSet(ts,StatsMonitor,&user,NULL);CHKERRQ(ierr);
   }
   if (output) {
-    ierr = TSMonitorSet(ts,OutputMonitor,&user,PETSC_NULL);CHKERRQ(ierr);
+    ierr = TSMonitorSet(ts,OutputMonitor,&user,NULL);CHKERRQ(ierr);
   }
 
   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
   ierr = VecDuplicate(U,&user.X0);CHKERRQ(ierr);
   ierr = VecCopy(U,user.X0);CHKERRQ(ierr);
 #if PETSC_VERSION_LE(3,3,0)
-  ierr = TSSolve(ts,U,PETSC_NULL);CHKERRQ(ierr);
+  ierr = TSSolve(ts,U,NULL);CHKERRQ(ierr);
 #else
   ierr = TSSolve(ts,U);CHKERRQ(ierr);
 #endif

demo/ClassicalShell.c

+/*
+  keywords: steady, vector, linear
+ */
 #include "petiga.h"
 
 typedef struct {

demo/ElasticRod.c

   user.E   = 1.0;
 
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","ElasticRod Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-density","Density",__FILE__,user.rho,&user.rho,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsReal("-Young_modulus","Young modulus",__FILE__,user.E,&user.E,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-density","Density",__FILE__,user.rho,&user.rho,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsReal("-Young_modulus","Young modulus",__FILE__,user.E,&user.E,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   IGA iga;

demo/Elasticity3D.c

+/*
+  This code solves the 3D elasticity equations given the Lame
+  constants and subject to Dirichlet boundary conditions.
+
+  keywords: steady, vector, linear
+ */
 #include "petiga.h"
 
 typedef struct {
   PetscReal lambda = user->lambda;
   PetscReal mu = user->mu;
 
-  PetscReal (*N1)[3] = (PetscReal (*)[3]) p->shape[1];
+  const PetscReal (*N1)[3];
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+
   PetscInt a,b,nen=p->nen;
   PetscScalar (*Kl)[3][nen][3] = (PetscScalar (*)[3][nen][3])K;
-
   for (a=0; a<nen; a++) {
     PetscReal Na_x = N1[a][0];
     PetscReal Na_y = N1[a][1];

demo/HyperElasticity.c

 #include "petscblaslapack.h"
 
 /*
-
   This code implements a HyperElastic material model in the context of
   large deformation elasticity. Implementation credit goes to students
   of the 2012 summer course `Nonlinear Finite Element Analysis' given
 #undef __FUNCT__
 #define __FUNCT__ "NeoHookeanModel"
 void NeoHookeanModel(IGAPoint pnt, const PetscScalar *U, PetscScalar (*F)[3], PetscScalar (*S)[3], PetscScalar (*D)[6], void *ctx)
-{    
+{
   AppCtx *user = (AppCtx *)ctx;
 
   PetscReal lambda = user->lambda;
   // Interpolate the solution and gradient given U
   PetscScalar u[3];
   PetscScalar grad_u[3][3];
-  IGAPointInterpolate(pnt,0,U,&u[0]);
-  IGAPointInterpolate(pnt,1,U,&grad_u[0][0]);
+  IGAPointFormValue(pnt,U,&u[0]);
+  IGAPointFormGrad (pnt,U,&grad_u[0][0]);
 
   // F = I + u_{i,j}
-  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2]; 
-  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2]; 
-  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2]; 
+  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2];
+  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2];
+  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2];
 
-  // Finv 
+  // Finv
   PetscScalar Finv[3][3],J,Jinv;
   J  = F[0][0]*(F[1][1]*F[2][2]-F[1][2]*F[2][1]);
   J -= F[0][1]*(F[1][0]*F[2][2]-F[1][2]*F[2][0]);
   J += F[0][2]*(F[1][0]*F[2][1]-F[1][1]*F[2][0]);
   Jinv = 1./J;
-  Finv[0][0] =  (F[1][1]*F[2][2]-F[2][1]*F[1][2])*Jinv; 
+  Finv[0][0] =  (F[1][1]*F[2][2]-F[2][1]*F[1][2])*Jinv;
   Finv[1][0] = -(F[1][0]*F[2][2]-F[1][2]*F[2][0])*Jinv;
   Finv[2][0] =  (F[1][0]*F[2][1]-F[2][0]*F[1][1])*Jinv;
   Finv[0][1] = -(F[0][1]*F[2][2]-F[0][2]*F[2][1])*Jinv;
-  Finv[1][1] =  (F[0][0]*F[2][2]-F[0][2]*F[2][0])*Jinv; 
+  Finv[1][1] =  (F[0][0]*F[2][2]-F[0][2]*F[2][0])*Jinv;
   Finv[2][1] = -(F[0][0]*F[2][1]-F[2][0]*F[0][1])*Jinv;
   Finv[0][2] =  (F[0][1]*F[1][2]-F[0][2]*F[1][1])*Jinv;
   Finv[1][2] = -(F[0][0]*F[1][2]-F[1][0]*F[0][2])*Jinv;
-  Finv[2][2] =  (F[0][0]*F[1][1]-F[1][0]*F[0][1])*Jinv; 
+  Finv[2][2] =  (F[0][0]*F[1][1]-F[1][0]*F[0][1])*Jinv;
 
   // C^-1 = (F^T F)^-1 = F^-1 F^-T
   PetscScalar Cinv[3][3];
   S[2][0] = temp*Cinv[2][0] + mu*(-Cinv[2][0]);
   S[2][1] = temp*Cinv[2][1] + mu*(-Cinv[2][1]);
   S[2][2] = temp*Cinv[2][2] + mu*(1.0-Cinv[2][2]);
-    
+
   // C_abcd=lambda*J^2*Cinv_ab*Cinv_cd+[2*miu-lambda(J^2-1)]*0.5(Cinv_ac*Cinv_bd+Cinv_ad*Cinv_bc)
   PetscScalar temp1=lambda*J*J;
   PetscScalar temp2=2*mu-lambda*(J*J-1);
   D[0][0]=temp1*Cinv[0][0]*Cinv[0][0]+temp2*0.5*(Cinv[0][0]*Cinv[0][0]+Cinv[0][0]*Cinv[0][0]);
-  D[0][1]=temp1*Cinv[0][0]*Cinv[1][1]+temp2*0.5*(Cinv[0][1]*Cinv[0][1]+Cinv[0][1]*Cinv[0][1]); 
-  D[0][2]=temp1*Cinv[0][0]*Cinv[2][2]+temp2*0.5*(Cinv[0][2]*Cinv[0][2]+Cinv[0][2]*Cinv[0][2]); 
-  D[0][3]=temp1*Cinv[0][0]*Cinv[0][1]+temp2*0.5*(Cinv[0][0]*Cinv[0][1]+Cinv[0][1]*Cinv[0][0]); 
-  D[0][4]=temp1*Cinv[0][0]*Cinv[1][2]+temp2*0.5*(Cinv[0][1]*Cinv[0][2]+Cinv[0][2]*Cinv[0][1]); 
-  D[0][5]=temp1*Cinv[0][0]*Cinv[0][2]+temp2*0.5*(Cinv[0][0]*Cinv[0][2]+Cinv[0][2]*Cinv[0][0]); 
-  D[1][1]=temp1*Cinv[1][1]*Cinv[1][1]+temp2*0.5*(Cinv[1][1]*Cinv[1][1]+Cinv[1][1]*Cinv[1][1]); 
-  D[1][2]=temp1*Cinv[1][1]*Cinv[2][2]+temp2*0.5*(Cinv[1][2]*Cinv[1][2]+Cinv[1][2]*Cinv[1][2]); 
-  D[1][3]=temp1*Cinv[1][1]*Cinv[0][1]+temp2*0.5*(Cinv[1][0]*Cinv[1][1]+Cinv[1][1]*Cinv[1][0]); 
-  D[1][4]=temp1*Cinv[1][1]*Cinv[1][2]+temp2*0.5*(Cinv[1][1]*Cinv[1][2]+Cinv[1][2]*Cinv[1][1]); 
-  D[1][5]=temp1*Cinv[1][1]*Cinv[0][2]+temp2*0.5*(Cinv[1][0]*Cinv[1][2]+Cinv[1][2]*Cinv[1][0]); 
-  D[2][2]=temp1*Cinv[2][2]*Cinv[2][2]+temp2*0.5*(Cinv[2][2]*Cinv[2][2]+Cinv[2][2]*Cinv[2][2]); 
-  D[2][3]=temp1*Cinv[2][2]*Cinv[0][1]+temp2*0.5*(Cinv[2][0]*Cinv[2][1]+Cinv[2][1]*Cinv[2][0]); 
-  D[2][4]=temp1*Cinv[2][2]*Cinv[1][2]+temp2*0.5*(Cinv[2][1]*Cinv[2][2]+Cinv[2][2]*Cinv[2][1]); 
-  D[2][5]=temp1*Cinv[2][2]*Cinv[0][2]+temp2*0.5*(Cinv[2][0]*Cinv[2][2]+Cinv[2][2]*Cinv[2][0]); 
-  D[3][3]=temp1*Cinv[0][1]*Cinv[0][1]+temp2*0.5*(Cinv[0][0]*Cinv[1][1]+Cinv[0][1]*Cinv[1][0]); 
-  D[3][4]=temp1*Cinv[0][1]*Cinv[1][2]+temp2*0.5*(Cinv[0][1]*Cinv[1][2]+Cinv[0][2]*Cinv[1][1]); 
-  D[3][5]=temp1*Cinv[0][1]*Cinv[0][2]+temp2*0.5*(Cinv[0][0]*Cinv[1][2]+Cinv[0][2]*Cinv[1][0]); 
-  D[4][4]=temp1*Cinv[1][2]*Cinv[1][2]+temp2*0.5*(Cinv[1][1]*Cinv[2][2]+Cinv[1][2]*Cinv[2][1]); 
-  D[4][5]=temp1*Cinv[1][2]*Cinv[0][2]+temp2*0.5*(Cinv[1][0]*Cinv[2][2]+Cinv[1][2]*Cinv[2][0]); 
-  D[5][5]=temp1*Cinv[0][2]*Cinv[0][2]+temp2*0.5*(Cinv[0][0]*Cinv[2][2]+Cinv[0][2]*Cinv[2][0]); 
+  D[0][1]=temp1*Cinv[0][0]*Cinv[1][1]+temp2*0.5*(Cinv[0][1]*Cinv[0][1]+Cinv[0][1]*Cinv[0][1]);
+  D[0][2]=temp1*Cinv[0][0]*Cinv[2][2]+temp2*0.5*(Cinv[0][2]*Cinv[0][2]+Cinv[0][2]*Cinv[0][2]);
+  D[0][3]=temp1*Cinv[0][0]*Cinv[0][1]+temp2*0.5*(Cinv[0][0]*Cinv[0][1]+Cinv[0][1]*Cinv[0][0]);
+  D[0][4]=temp1*Cinv[0][0]*Cinv[1][2]+temp2*0.5*(Cinv[0][1]*Cinv[0][2]+Cinv[0][2]*Cinv[0][1]);
+  D[0][5]=temp1*Cinv[0][0]*Cinv[0][2]+temp2*0.5*(Cinv[0][0]*Cinv[0][2]+Cinv[0][2]*Cinv[0][0]);
+  D[1][1]=temp1*Cinv[1][1]*Cinv[1][1]+temp2*0.5*(Cinv[1][1]*Cinv[1][1]+Cinv[1][1]*Cinv[1][1]);
+  D[1][2]=temp1*Cinv[1][1]*Cinv[2][2]+temp2*0.5*(Cinv[1][2]*Cinv[1][2]+Cinv[1][2]*Cinv[1][2]);
+  D[1][3]=temp1*Cinv[1][1]*Cinv[0][1]+temp2*0.5*(Cinv[1][0]*Cinv[1][1]+Cinv[1][1]*Cinv[1][0]);
+  D[1][4]=temp1*Cinv[1][1]*Cinv[1][2]+temp2*0.5*(Cinv[1][1]*Cinv[1][2]+Cinv[1][2]*Cinv[1][1]);
+  D[1][5]=temp1*Cinv[1][1]*Cinv[0][2]+temp2*0.5*(Cinv[1][0]*Cinv[1][2]+Cinv[1][2]*Cinv[1][0]);
+  D[2][2]=temp1*Cinv[2][2]*Cinv[2][2]+temp2*0.5*(Cinv[2][2]*Cinv[2][2]+Cinv[2][2]*Cinv[2][2]);
+  D[2][3]=temp1*Cinv[2][2]*Cinv[0][1]+temp2*0.5*(Cinv[2][0]*Cinv[2][1]+Cinv[2][1]*Cinv[2][0]);
+  D[2][4]=temp1*Cinv[2][2]*Cinv[1][2]+temp2*0.5*(Cinv[2][1]*Cinv[2][2]+Cinv[2][2]*Cinv[2][1]);
+  D[2][5]=temp1*Cinv[2][2]*Cinv[0][2]+temp2*0.5*(Cinv[2][0]*Cinv[2][2]+Cinv[2][2]*Cinv[2][0]);
+  D[3][3]=temp1*Cinv[0][1]*Cinv[0][1]+temp2*0.5*(Cinv[0][0]*Cinv[1][1]+Cinv[0][1]*Cinv[1][0]);
+  D[3][4]=temp1*Cinv[0][1]*Cinv[1][2]+temp2*0.5*(Cinv[0][1]*Cinv[1][2]+Cinv[0][2]*Cinv[1][1]);
+  D[3][5]=temp1*Cinv[0][1]*Cinv[0][2]+temp2*0.5*(Cinv[0][0]*Cinv[1][2]+Cinv[0][2]*Cinv[1][0]);
+  D[4][4]=temp1*Cinv[1][2]*Cinv[1][2]+temp2*0.5*(Cinv[1][1]*Cinv[2][2]+Cinv[1][2]*Cinv[2][1]);
+  D[4][5]=temp1*Cinv[1][2]*Cinv[0][2]+temp2*0.5*(Cinv[1][0]*Cinv[2][2]+Cinv[1][2]*Cinv[2][0]);
+  D[5][5]=temp1*Cinv[0][2]*Cinv[0][2]+temp2*0.5*(Cinv[0][0]*Cinv[2][2]+Cinv[0][2]*Cinv[2][0]);
   D[1][0]=D[0][1];
   D[2][0]=D[0][2];
   D[3][0]=D[0][3];
   D[5][2]=D[2][5];
   D[4][3]=D[3][4];
   D[5][3]=D[3][5];
-  D[5][4]=D[4][5];  
+  D[5][4]=D[4][5];
 }
 
 #undef __FUNCT__
 #define __FUNCT__ "StVenantModel"
 void StVenantModel(IGAPoint pnt, const PetscScalar *U, PetscScalar (*F)[3], PetscScalar (*S)[3], PetscScalar (*D)[6], void *ctx)
-{    
+{
   AppCtx *user = (AppCtx *)ctx;
 
   PetscReal lambda = user->lambda;
   // Interpolate the solution and gradient given U
   PetscScalar u[3];
   PetscScalar grad_u[3][3];
-  IGAPointInterpolate(pnt,0,U,&u[0]);
-  IGAPointInterpolate(pnt,1,U,&grad_u[0][0]);
+  IGAPointFormValue(pnt,U,&u[0]);
+  IGAPointFormGrad (pnt,U,&grad_u[0][0]);
 
   // F = I + u_{i,j}
-  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2]; 
-  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2]; 
-  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2]; 
+  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2];
+  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2];
+  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2];
 
   // C^-1 = (F^T F)^-1 = F^-1 F^-T
   PetscScalar C[3][3];
   S[2][0] = mu*C[2][0];
   S[2][1] = mu*C[2][1];
   S[2][2] = temp + mu*(C[2][2]-1);
-    
-  //C_abcd = lambda*delta_ab*delta_cd+2*mu*0.5*(delta_ac*delta_bd+delta_ad*delta_bc) 
+
+  //C_abcd = lambda*delta_ab*delta_cd+2*mu*0.5*(delta_ac*delta_bd+delta_ad*delta_bc)
   //Elasticity material tensor
   D[0][0]=lambda+2*mu;
   D[0][1]=lambda;
   D[5][2]=D[2][5];
   D[4][3]=D[3][4];
   D[5][3]=D[3][5];
-  D[5][4]=D[4][5];  
+  D[5][4]=D[4][5];
 }
 
 #undef __FUNCT__
 #define __FUNCT__ "MooneyRivlinModel1"
 void MooneyRivlinModel1(IGAPoint pnt, const PetscScalar *U, PetscScalar (*F)[3], PetscScalar (*S)[3], PetscScalar (*D)[6], void *ctx)
-{    
+{
   AppCtx *user = (AppCtx *)ctx;
 
   PetscReal a = user->a;
   // Interpolate the solution and gradient given U
   PetscScalar u[3];
   PetscScalar grad_u[3][3];
-  IGAPointInterpolate(pnt,0,U,&u[0]);
-  IGAPointInterpolate(pnt,1,U,&grad_u[0][0]);
+  IGAPointFormValue(pnt,U,&u[0]);
+  IGAPointFormGrad (pnt,U,&grad_u[0][0]);
 
   // F = I + u_{i,j}
-  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2]; 
-  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2]; 
-  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2]; 
+  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2];
+  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2];
+  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2];
 
-  // Finv 
+  // Finv
   PetscScalar Finv[3][3],J,Jinv;
   J  = F[0][0]*(F[1][1]*F[2][2]-F[1][2]*F[2][1]);
   J -= F[0][1]*(F[1][0]*F[2][2]-F[1][2]*F[2][0]);
   J += F[0][2]*(F[1][0]*F[2][1]-F[1][1]*F[2][0]);
   Jinv = 1./J;
-  Finv[0][0] =  (F[1][1]*F[2][2]-F[2][1]*F[1][2])*Jinv; 
+  Finv[0][0] =  (F[1][1]*F[2][2]-F[2][1]*F[1][2])*Jinv;
   Finv[1][0] = -(F[1][0]*F[2][2]-F[1][2]*F[2][0])*Jinv;
   Finv[2][0] =  (F[1][0]*F[2][1]-F[2][0]*F[1][1])*Jinv;
   Finv[0][1] = -(F[0][1]*F[2][2]-F[0][2]*F[2][1])*Jinv;
-  Finv[1][1] =  (F[0][0]*F[2][2]-F[0][2]*F[2][0])*Jinv; 
+  Finv[1][1] =  (F[0][0]*F[2][2]-F[0][2]*F[2][0])*Jinv;
   Finv[2][1] = -(F[0][0]*F[2][1]-F[2][0]*F[0][1])*Jinv;
   Finv[0][2] =  (F[0][1]*F[1][2]-F[0][2]*F[1][1])*Jinv;
   Finv[1][2] = -(F[0][0]*F[1][2]-F[1][0]*F[0][2])*Jinv;
-  Finv[2][2] =  (F[0][0]*F[1][1]-F[1][0]*F[0][1])*Jinv; 
+  Finv[2][2] =  (F[0][0]*F[1][1]-F[1][0]*F[0][1])*Jinv;
 
   // C^-1 = (F^T F)^-1 = F^-1 F^-T
   PetscScalar Cinv[3][3];
   S[2][0]=(J*J*(2*c+b*temp)-d)*Cinv[2][0]-b*J*J*Cs[2][0];
   S[2][1]=(J*J*(2*c+b*temp)-d)*Cinv[2][1]-b*J*J*Cs[2][1];
   S[2][2]=2*a+(J*J*(2*c+b*temp)-d)*Cinv[2][2]-b*J*J*Cs[2][2];
-    
+
   //Dijkl = J^2(4c+2b*tr(Cinv))*Cinvij*Cinvkl-2bJ^2(Cinv^2ij*Cinvkl+Cinvij*Cinv^2kl-0.5*(Cinv^2li*Cinvjk+Cinvli*Cinv^2jk+Cinv^2lj*Cinvik+Cinvlj*Cinv^2ik))-(J^2*(4c+2btr(Cinv))-2d)*0.5*(Cinvik*Cinvlj+Cinvli*Cinvjk)
   PetscScalar temp1= J*J*(4*c+2*b*temp);
   PetscScalar temp2=2*b*J*J;
 #undef __FUNCT__
 #define __FUNCT__ "MooneyRivlinModel2"
 void MooneyRivlinModel2(IGAPoint pnt, const PetscScalar *U, PetscScalar (*F)[3], PetscScalar (*S)[3], PetscScalar (*D)[6], void *ctx)
-{    
+{
   AppCtx *user = (AppCtx *)ctx;
 
   PetscReal c1 = user->c1;
   // Interpolate the solution and gradient given U
   PetscScalar u[3];
   PetscScalar grad_u[3][3];
-  IGAPointInterpolate(pnt,0,U,&u[0]);
-  IGAPointInterpolate(pnt,1,U,&grad_u[0][0]);
+  IGAPointFormValue(pnt,U,&u[0]);
+  IGAPointFormGrad (pnt,U,&grad_u[0][0]);
 
   // F = I + u_{i,j}
-  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2]; 
-  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2]; 
-  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2]; 
+  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2];
+  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2];
+  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2];
 
-  // Finv 
+  // Finv
   PetscScalar Finv[3][3],J,Jinv;
   J  = F[0][0]*(F[1][1]*F[2][2]-F[1][2]*F[2][1]);
   J -= F[0][1]*(F[1][0]*F[2][2]-F[1][2]*F[2][0]);
   J += F[0][2]*(F[1][0]*F[2][1]-F[1][1]*F[2][0]);
   Jinv = 1./J;
-  Finv[0][0] =  (F[1][1]*F[2][2]-F[2][1]*F[1][2])*Jinv; 
+  Finv[0][0] =  (F[1][1]*F[2][2]-F[2][1]*F[1][2])*Jinv;
   Finv[1][0] = -(F[1][0]*F[2][2]-F[1][2]*F[2][0])*Jinv;
   Finv[2][0] =  (F[1][0]*F[2][1]-F[2][0]*F[1][1])*Jinv;
   Finv[0][1] = -(F[0][1]*F[2][2]-F[0][2]*F[2][1])*Jinv;
-  Finv[1][1] =  (F[0][0]*F[2][2]-F[0][2]*F[2][0])*Jinv; 
+  Finv[1][1] =  (F[0][0]*F[2][2]-F[0][2]*F[2][0])*Jinv;
   Finv[2][1] = -(F[0][0]*F[2][1]-F[2][0]*F[0][1])*Jinv;
   Finv[0][2] =  (F[0][1]*F[1][2]-F[0][2]*F[1][1])*Jinv;
   Finv[1][2] = -(F[0][0]*F[1][2]-F[1][0]*F[0][2])*Jinv;
-  Finv[2][2] =  (F[0][0]*F[1][1]-F[1][0]*F[0][1])*Jinv; 
+  Finv[2][2] =  (F[0][0]*F[1][1]-F[1][0]*F[0][1])*Jinv;
 
   // C^-1 = (F^T F)^-1 = F^-1 F^-T
   PetscScalar Cinv[3][3];
   trCs += C[2][0]*C[0][2] + C[2][1]*C[1][2] + C[2][2]*C[2][2];
 
   PetscScalar sinva = 0.5*(trC*trC - trCs);
-  
+
   PetscScalar temp = pow(J,-2/3);
   PetscScalar tempD = (c1*trC+2*c2*temp*sinva);
-	
+
   //Sij=2*a*deltaij+(J^2*(2*c+b*tr(Cinv))-d)Cinvij-b*J^2*Csij
   PetscScalar Siso[3][3];
   Siso[0][0]=2*temp*(c1+c2*temp*(trC-C[0][0])-(1/3)*tempD*Cinv[0][0]);
   Dvol[4][4]=temp1*Cinv[1][2]*Cinv[1][2]-temp2*0.5*(Cinv[1][1]*Cinv[2][2]+Cinv[2][1]*Cinv[2][1]);
   Dvol[4][5]=temp1*Cinv[1][2]*Cinv[0][2]-temp2*0.5*(Cinv[1][0]*Cinv[2][2]+Cinv[2][1]*Cinv[2][0]);
   Dvol[5][5]=temp1*Cinv[0][2]*Cinv[0][2]-temp2*0.5*(Cinv[0][0]*Cinv[2][2]+Cinv[2][0]*Cinv[2][0]);
-  
+
 
   //Dijkl = J^2(4c+2b*tr(Cinv))*Cinvij*Cinvkl-2bJ^2(Cinv^2ij*Cinvkl+Cinvij*Cinv^2kl-0.5*(Cinv^2li*Cinvjk+Cinvli*Cinv^2jk+Cinv^2lj*Cinvik+Cinvlj*Cinv^2ik))-(J^2*(4c+2btr(Cinv))-2d)*0.5*(Cinvik*Cinvlj+Cinvli*Cinvjk)
   D[0][0]=Diso[0][0]+Dvol[0][0];
 #undef __FUNCT__
 #define __FUNCT__ "GeneralModel"
 void GeneralModel(IGAPoint pnt, const PetscScalar *U, PetscScalar (*F)[3], PetscScalar (*S)[3], PetscScalar (*D)[6], void *ctx)
-{    
+{
   AppCtx *user = (AppCtx *)ctx;
 
   PetscReal c1 = user->c1;
   // Interpolate the solution and gradient given U
   PetscScalar u[3];
   PetscScalar grad_u[3][3];
-  IGAPointInterpolate(pnt,0,U,&u[0]);
-  IGAPointInterpolate(pnt,1,U,&grad_u[0][0]);
+  IGAPointFormValue(pnt,U,&u[0]);
+  IGAPointFormGrad (pnt,U,&grad_u[0][0]);
 
   // F = I + u_{i,j}
-  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2]; 
-  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2]; 
-  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2]; 
+  F[0][0] = 1+grad_u[0][0]; F[0][1] =   grad_u[0][1];  F[0][2] =   grad_u[0][2];
+  F[1][0] =   grad_u[1][0]; F[1][1] = 1+grad_u[1][1];  F[1][2] =   grad_u[1][2];
+  F[2][0] =   grad_u[2][0]; F[2][1] =   grad_u[2][1];  F[2][2] = 1+grad_u[2][2];
 
-  // Finv 
+  // Finv
   PetscScalar Finv[3][3],J,Jinv;
   J  = F[0][0]*(F[1][1]*F[2][2]-F[1][2]*F[2][1]);
   J -= F[0][1]*(F[1][0]*F[2][2]-F[1][2]*F[2][0]);
   J += F[0][2]*(F[1][0]*F[2][1]-F[1][1]*F[2][0]);
   Jinv = 1./J;
-  Finv[0][0] =  (F[1][1]*F[2][2]-F[2][1]*F[1][2])*Jinv; 
+  Finv[0][0] =  (F[1][1]*F[2][2]-F[2][1]*F[1][2])*Jinv;
   Finv[1][0] = -(F[1][0]*F[2][2]-F[1][2]*F[2][0])*Jinv;
   Finv[2][0] =  (F[1][0]*F[2][1]-F[2][0]*F[1][1])*Jinv;
   Finv[0][1] = -(F[0][1]*F[2][2]-F[0][2]*F[2][1])*Jinv;
-  Finv[1][1] =  (F[0][0]*F[2][2]-F[0][2]*F[2][0])*Jinv; 
+  Finv[1][1] =  (F[0][0]*F[2][2]-F[0][2]*F[2][0])*Jinv;
   Finv[2][1] = -(F[0][0]*F[2][1]-F[2][0]*F[0][1])*Jinv;
   Finv[0][2] =  (F[0][1]*F[1][2]-F[0][2]*F[1][1])*Jinv;
   Finv[1][2] = -(F[0][0]*F[1][2]-F[1][0]*F[0][2])*Jinv;
   P[5][5] = 2*tempJ*(0.5-(1/3)*Cinv[0][2]*C[0][2]);
 
   PetscScalar trC = C[0][0]+C[1][1]+C[2][2];
-  
+
   PetscScalar Shat[6];
   Shat[0] = 2*(c1+c2*tempJ*(trC-C[0][0])); //Shat[0][0]
   Shat[1] = 2*(c1+c2*tempJ*(trC-C[1][1])); //Shat[1][1]
   Siso[3]=P[3][0]*Shat[0] + P[3][1]*Shat[1] + P[3][2]*Shat[2] + P[3][3]*Shat[3] + P[3][4]*Shat[4] + P[3][5]*Shat[5];
   Siso[4]=P[4][0]*Shat[0] + P[4][1]*Shat[1] + P[4][2]*Shat[2] + P[4][3]*Shat[3] + P[4][4]*Shat[4] + P[4][5]*Shat[5];
   Siso[5]=P[5][0]*Shat[0] + P[5][1]*Shat[1] + P[5][2]*Shat[2] + P[5][3]*Shat[3] + P[5][4]*Shat[4] + P[5][5]*Shat[5];
-	
+
   PetscScalar tempdJ = kappa*(J-1);
   PetscScalar tempd2J2 = kappa;
   PetscScalar Svol[6];
   Dhat[4][3]=Dhat[3][4];
   Dhat[5][3]=Dhat[3][5];
   Dhat[5][4]=Dhat[4][5];
-  
+
 
   PetscScalar Diso[6][6];
-  Diso[0][0] = P[0][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[0][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[0][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[0][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[0][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[0][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]); 
+  Diso[0][0] = P[0][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[0][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[0][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[0][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[0][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[0][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]);
   Diso[0][1] = P[1][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[1][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[1][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[1][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[1][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[1][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]);
   Diso[0][2] = P[2][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[2][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[2][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[2][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[2][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[2][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]);
-  Diso[0][3] = P[3][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[3][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[3][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[3][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[3][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[3][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]); 
-  Diso[0][4] = P[4][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[4][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[4][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[4][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[4][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[4][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]); 
+  Diso[0][3] = P[3][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[3][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[3][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[3][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[3][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[3][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]);
+  Diso[0][4] = P[4][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[4][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[4][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[4][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[4][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[4][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]);
   Diso[0][5] = P[5][0]*(Dhat[0][0]*P[0][0] + Dhat[1][0]*P[0][1] + Dhat[2][0]*P[0][2] + Dhat[3][0]*P[0][3] + Dhat[4][0]*P[0][4] + Dhat[5][0]*P[0][5]) + P[5][1]*(Dhat[0][1]*P[0][0] + Dhat[1][1]*P[0][1] + Dhat[2][1]*P[0][2] + Dhat[3][1]*P[0][3] + Dhat[4][1]*P[0][4] + Dhat[5][1]*P[0][5]) + P[5][2]*(Dhat[0][2]*P[0][0] + Dhat[1][2]*P[0][1] + Dhat[2][2]*P[0][2] + Dhat[3][2]*P[0][3] + Dhat[4][2]*P[0][4] + Dhat[5][2]*P[0][5]) + P[5][3]*(Dhat[0][3]*P[0][0] + Dhat[1][3]*P[0][1] + Dhat[2][3]*P[0][2] + Dhat[3][3]*P[0][3] + Dhat[4][3]*P[0][4] + Dhat[5][3]*P[0][5]) + P[5][4]*(Dhat[0][4]*P[0][0] + Dhat[1][4]*P[0][1] + Dhat[2][4]*P[0][2] + Dhat[3][4]*P[0][3] + Dhat[4][4]*P[0][4] + Dhat[5][4]*P[0][5]) + P[5][5]*(Dhat[0][5]*P[0][0] + Dhat[1][5]*P[0][1] + Dhat[2][5]*P[0][2] + Dhat[3][5]*P[0][3] + Dhat[4][5]*P[0][4] + Dhat[5][5]*P[0][5]);
-  Diso[1][1] = P[1][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[1][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[1][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[1][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[1][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[1][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]); 
-  Diso[1][2] = P[2][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[2][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[2][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[2][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[2][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[2][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]); 
-  Diso[1][3] = P[3][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[3][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[3][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[3][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[3][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[3][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]); 
-  Diso[1][4] = P[4][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[4][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[4][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[4][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[4][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[4][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]); 
+  Diso[1][1] = P[1][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[1][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[1][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[1][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[1][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[1][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]);
+  Diso[1][2] = P[2][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[2][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[2][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[2][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[2][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[2][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]);
+  Diso[1][3] = P[3][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[3][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[3][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[3][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[3][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[3][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]);
+  Diso[1][4] = P[4][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[4][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[4][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[4][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[4][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[4][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]);
   Diso[1][5] = P[5][0]*(Dhat[0][0]*P[1][0] + Dhat[1][0]*P[1][1] + Dhat[2][0]*P[1][2] + Dhat[3][0]*P[1][3] + Dhat[4][0]*P[1][4] + Dhat[5][0]*P[1][5]) + P[5][1]*(Dhat[0][1]*P[1][0] + Dhat[1][1]*P[1][1] + Dhat[2][1]*P[1][2] + Dhat[3][1]*P[1][3] + Dhat[4][1]*P[1][4] + Dhat[5][1]*P[1][5]) + P[5][2]*(Dhat[0][2]*P[1][0] + Dhat[1][2]*P[1][1] + Dhat[2][2]*P[1][2] + Dhat[3][2]*P[1][3] + Dhat[4][2]*P[1][4] + Dhat[5][2]*P[1][5]) + P[5][3]*(Dhat[0][3]*P[1][0] + Dhat[1][3]*P[1][1] + Dhat[2][3]*P[1][2] + Dhat[3][3]*P[1][3] + Dhat[4][3]*P[1][4] + Dhat[5][3]*P[1][5]) + P[5][4]*(Dhat[0][4]*P[1][0] + Dhat[1][4]*P[1][1] + Dhat[2][4]*P[1][2] + Dhat[3][4]*P[1][3] + Dhat[4][4]*P[1][4] + Dhat[5][4]*P[1][5]) + P[5][5]*(Dhat[0][5]*P[1][0] + Dhat[1][5]*P[1][1] + Dhat[2][5]*P[1][2] + Dhat[3][5]*P[1][3] + Dhat[4][5]*P[1][4] + Dhat[5][5]*P[1][5]);
-  Diso[2][2] = P[2][0]*(Dhat[0][0]*P[2][0] + Dhat[1][0]*P[2][1] + Dhat[2][0]*P[2][2] + Dhat[3][0]*P[2][3] + Dhat[4][0]*P[2][4] + Dhat[5][0]*P[2][5]) + P[2][1]*(Dhat[0][1]*P[2][0] + Dhat[1][1]*P[2][1] + Dhat[2][1]*P[2][2] + Dhat[3][1]*P[2][3] + Dhat[4][1]*P[2][4] + Dhat[5][1]*P[2][5]) + P[2][2]*(Dhat[0][2]*P[2][0] + Dhat[1][2]*P[2][1] + Dhat[2][2]*P[2][2] + Dhat[3][2]*P[2][3] + Dhat[4][2]*P[2][4] + Dhat[5][2]*P[2][5]) + P[2][3]*(Dhat[0][3]*P[2][0] + Dhat[1][3]*P[2][1] + Dhat[2][3]*P[2][2] + Dhat[3][3]*P[2][3] + Dhat[4][3]*P[2][4] + Dhat[5][3]*P[2][5]) + P[2][4]*(Dhat[0][4]*P[2][0] + Dhat[1][4]*P[2][1] + Dhat[2][4]*P[2][2] + Dhat[3][4]*P[2][3] + Dhat[4][4]*P[2][4] + Dhat[5][4]*P[2][5]) + P[2][5]*(Dhat[0][5]*P[2][0] + Dhat[1][5]*P[2][1] + Dhat[2][5]*P[2][2] + Dhat[3][5]*P[2][3] + Dhat[4][5]*P[2][4] + Dhat[5][5]*P[2][5]); 
-  Diso[2][3] = P[3][0]*(Dhat[0][0]*P[2][0] + Dhat[1][0]*P[2][1] + Dhat[2][0]*P[2][2] + Dhat[3][0]*P[2][3] + Dhat[4][0]*P[2][4] + Dhat[5][0]*P[2][5]) + P[3][1]*(Dhat[0][1]*P[2][0] + Dhat[1][1]*P[2][1] + Dhat[2][1]*P[2][2] + Dhat[3][1]*P[2][3] + Dhat[4][1]*P[2][4] + Dhat[5][1]*P[2][5]) + P[3][2]*(Dhat[0][2]*P[2][0] + Dhat[1][2]*P[2][1] + Dhat[2][2]*P[2][2] + Dhat[3][2]*P[2][3] + Dhat[4][2]*P[2][4] + Dhat[5][2]*P[2][5]) + P[3][3]*(Dhat[0][3]*P[2][0] + Dhat[1][3]*P[2][1] + Dhat[2][3]*P[2][2] + Dhat[3][3]*P[2][3] + Dhat[4][3]*P[2][4] + Dhat[5][3]*P[2][5]) + P[3][4]*(Dhat[0][4]*P[2][0] + Dhat[1][4]*P[2][1] + Dhat[2][4]*P[2][2] + Dhat[3][4]*P[2][3] + Dhat[4][4]*P[2][4] + Dhat[5][4]*P[2][5]) + P[3][5]*(Dhat[0][5]*P[2][0] + Dhat[1][5]*P[2][1] + Dhat[2][5]*P[2][2] + Dhat[3][5]*P[2][3] + Dhat[4][5]*P[2][4] + Dhat[5][5]*P[2][5]); 
-  Diso[2][4] = P[4][0]*(Dhat[0][0]*P[2][0] + Dhat[1][0]*P[2][1] + Dhat[2][0]*P[2][2] + Dhat[3][0]*P[2][3] + Dhat[4][0]*P[2][4] + Dhat[5][0]*P[2][5]) + P[4][1]*(Dhat[0][1]*P[2][0] + Dhat[1][1]*P[2][1] + Dhat[2][1]*P[2][2] + Dhat[3][1]*P[2][3] + Dhat[4][1]*P[2][4] + Dhat[5][1]*P[2][5]) + P[4][2]*(Dhat[0][2]*P[2][0] + Dhat[1][2]*P[2][1] + Dhat[2][2]*P[2][2] + Dhat[3][2]*P[2][3] + Dhat[4][2]*P[2][4] + Dhat[5][2]*P[2][5]) + P[4][3]*(Dhat[0][3]*P[2][0] + Dhat[1][3]*P[2][1] + Dhat[2][3]*P[2][2] + Dhat[3][3]*P[2][3] + Dhat[4][3]*P[2][4] + Dhat[5][3]*P[2][5]) + P[4][4]*(Dhat[0][4]*P[2][0] + Dhat[1][4]*P[2][1] + Dhat[2][4]*P[2][2] + Dhat[3][4]*P[2][3] + Dhat[4][4]*P[2][4] + Dhat[5][4]*P[2][5]) + P[4][5]*(Dhat[0][5]*P[2][0] + Dhat[1][5]*P[2][1] + Dhat[2][5]*P[2][2] + Dhat[3][5]*P[2][3] + Dhat[4][5]*P[2][4] + Dhat[5][5]*P[2][5]); 
+  Diso[2][2] = P[2][0]*(Dhat[0][0]*P[2][0] + Dhat[1][0]*P[2][1] + Dhat[2][0]*P[2][2] + Dhat[3][0]*P[2][3] + Dhat[4][0]*P[2][4] + Dhat[5][0]*P[2][5]) + P[2][1]*(Dhat[0][1]*P[2][0] + Dhat[1][1]*P[2][1] + Dhat[2][1]*P[2][2] + Dhat[3][1]*P[2][3] + Dhat[4][1]*P[2][4] + Dhat[5][1]*P[2][5]) + P[2][2]*(Dhat[0][2]*P[2][0] + Dhat[1][2]*P[2][1] + Dhat[2][2]*P[2][2] + Dhat[3][2]*P[2][3] + Dhat[4][2]*P[2][4] + Dhat[5][2]*P[2][5]) + P[2][3]*(Dhat[0][3]*P[2][0] + Dhat[1][3]*P[2][1] + Dhat[2][3]*P[2][2] + Dhat[3][3]*P[2][3] + Dhat[4][3]*P[2][4] + Dhat[5][3]*P[2][5]) + P[2][4]*(Dhat[0][4]*P[2][0] + Dhat[1][4]*P[2][1] + Dhat[2][4]*P[2][2] + Dhat[3][4]*P[2][3] + Dhat[4][4]*P[2][4] + Dhat[5][4]*P[2][5]) + P[2][5]*(Dhat[0][5]*P[2][0] + Dhat[1][5]*P[2][1] + Dhat[2][5]*P[2][2] + Dhat[3][5]*P[2][3] + Dhat[4][5]*P[2][4] + Dhat[5][5]*P[2][5]);
+  Diso[2][3] = P[3][0]*(Dhat[0][0]*P[2][0] + Dhat[1][0]*P[2][1] + Dhat[2][0]*P[2][2] + Dhat[3][0]*P[2][3] + Dhat[4][0]*P[2][4] + Dhat[5][0]*P[2][5]) + P[3][1]*(Dhat[0][1]*P[2][0] + Dhat[1][1]*P[2][1] + Dhat[2][1]*P[2][2] + Dhat[3][1]*P[2][3] + Dhat[4][1]*P[2][4] + Dhat[5][1]*P[2][5]) + P[3][2]*(Dhat[0][2]*P[2][0] + Dhat[1][2]*P[2][1] + Dhat[2][2]*P[2][2] + Dhat[3][2]*P[2][3] + Dhat[4][2]*P[2][4] + Dhat[5][2]*P[2][5]) + P[3][3]*(Dhat[0][3]*P[2][0] + Dhat[1][3]*P[2][1] + Dhat[2][3]*P[2][2] + Dhat[3][3]*P[2][3] + Dhat[4][3]*P[2][4] + Dhat[5][3]*P[2][5]) + P[3][4]*(Dhat[0][4]*P[2][0] + Dhat[1][4]*P[2][1] + Dhat[2][4]*P[2][2] + Dhat[3][4]*P[2][3] + Dhat[4][4]*P[2][4] + Dhat[5][4]*P[2][5]) + P[3][5]*(Dhat[0][5]*P[2][0] + Dhat[1][5]*P[2][1] + Dhat[2][5]*P[2][2] + Dhat[3][5]*P[2][3] + Dhat[4][5]*P[2][4] + Dhat[5][5]*P[2][5]);
+  Diso[2][4] = P[4][0]*(Dhat[0][0]*P[2][0] + Dhat[1][0]*P[2][1] + Dhat[2][0]*P[2][2] + Dhat[3][0]*P[2][3] + Dhat[4][0]*P[2][4] + Dhat[5][0]*P[2][5]) + P[4][1]*(Dhat[0][1]*P[2][0] + Dhat[1][1]*P[2][1] + Dhat[2][1]*P[2][2] + Dhat[3][1]*P[2][3] + Dhat[4][1]*P[2][4] + Dhat[5][1]*P[2][5]) + P[4][2]*(Dhat[0][2]*P[2][0] + Dhat[1][2]*P[2][1] + Dhat[2][2]*P[2][2] + Dhat[3][2]*P[2][3] + Dhat[4][2]*P[2][4] + Dhat[5][2]*P[2][5]) + P[4][3]*(Dhat[0][3]*P[2][0] + Dhat[1][3]*P[2][1] + Dhat[2][3]*P[2][2] + Dhat[3][3]*P[2][3] + Dhat[4][3]*P[2][4] + Dhat[5][3]*P[2][5]) + P[4][4]*(Dhat[0][4]*P[2][0] + Dhat[1][4]*P[2][1] + Dhat[2][4]*P[2][2] + Dhat[3][4]*P[2][3] + Dhat[4][4]*P[2][4] + Dhat[5][4]*P[2][5]) + P[4][5]*(Dhat[0][5]*P[2][0] + Dhat[1][5]*P[2][1] + Dhat[2][5]*P[2][2] + Dhat[3][5]*P[2][3] + Dhat[4][5]*P[2][4] + Dhat[5][5]*P[2][5]);
   Diso[2][5] = P[5][0]*(Dhat[0][0]*P[2][0] + Dhat[1][0]*P[2][1] + Dhat[2][0]*P[2][2] + Dhat[3][0]*P[2][3] + Dhat[4][0]*P[2][4] + Dhat[5][0]*P[2][5]) + P[5][1]*(Dhat[0][1]*P[2][0] + Dhat[1][1]*P[2][1] + Dhat[2][1]*P[2][2] + Dhat[3][1]*P[2][3] + Dhat[4][1]*P[2][4] + Dhat[5][1]*P[2][5]) + P[5][2]*(Dhat[0][2]*P[2][0] + Dhat[1][2]*P[2][1] + Dhat[2][2]*P[2][2] + Dhat[3][2]*P[2][3] + Dhat[4][2]*P[2][4] + Dhat[5][2]*P[2][5]) + P[5][3]*(Dhat[0][3]*P[2][0] + Dhat[1][3]*P[2][1] + Dhat[2][3]*P[2][2] + Dhat[3][3]*P[2][3] + Dhat[4][3]*P[2][4] + Dhat[5][3]*P[2][5]) + P[5][4]*(Dhat[0][4]*P[2][0] + Dhat[1][4]*P[2][1] + Dhat[2][4]*P[2][2] + Dhat[3][4]*P[2][3] + Dhat[4][4]*P[2][4] + Dhat[5][4]*P[2][5]) + P[5][5]*(Dhat[0][5]*P[2][0] + Dhat[1][5]*P[2][1] + Dhat[2][5]*P[2][2] + Dhat[3][5]*P[2][3] + Dhat[4][5]*P[2][4] + Dhat[5][5]*P[2][5]);
-  Diso[3][3] = P[3][0]*(Dhat[0][0]*P[3][0] + Dhat[1][0]*P[3][1] + Dhat[2][0]*P[3][2] + Dhat[3][0]*P[3][3] + Dhat[4][0]*P[3][4] + Dhat[5][0]*P[3][5]) + P[3][1]*(Dhat[0][1]*P[3][0] + Dhat[1][1]*P[3][1] + Dhat[2][1]*P[3][2] + Dhat[3][1]*P[3][3] + Dhat[4][1]*P[3][4] + Dhat[5][1]*P[3][5]) + P[3][2]*(Dhat[0][2]*P[3][0] + Dhat[1][2]*P[3][1] + Dhat[2][2]*P[3][2] + Dhat[3][2]*P[3][3] + Dhat[4][2]*P[3][4] + Dhat[5][2]*P[3][5]) + P[3][3]*(Dhat[0][3]*P[3][0] + Dhat[1][3]*P[3][1] + Dhat[2][3]*P[3][2] + Dhat[3][3]*P[3][3] + Dhat[4][3]*P[3][4] + Dhat[5][3]*P[3][5]) + P[3][4]*(Dhat[0][4]*P[3][0] + Dhat[1][4]*P[3][1] + Dhat[2][4]*P[3][2] + Dhat[3][4]*P[3][3] + Dhat[4][4]*P[3][4] + Dhat[5][4]*P[3][5]) + P[3][5]*(Dhat[0][5]*P[3][0] + Dhat[1][5]*P[3][1] + Dhat[2][5]*P[3][2] + Dhat[3][5]*P[3][3] + Dhat[4][5]*P[3][4] + Dhat[5][5]*P[3][5]); 
-  Diso[3][4] = P[4][0]*(Dhat[0][0]*P[3][0] + Dhat[1][0]*P[3][1] + Dhat[2][0]*P[3][2] + Dhat[3][0]*P[3][3] + Dhat[4][0]*P[3][4] + Dhat[5][0]*P[3][5]) + P[4][1]*(Dhat[0][1]*P[3][0] + Dhat[1][1]*P[3][1] + Dhat[2][1]*P[3][2] + Dhat[3][1]*P[3][3] + Dhat[4][1]*P[3][4] + Dhat[5][1]*P[3][5]) + P[4][2]*(Dhat[0][2]*P[3][0] + Dhat[1][2]*P[3][1] + Dhat[2][2]*P[3][2] + Dhat[3][2]*P[3][3] + Dhat[4][2]*P[3][4] + Dhat[5][2]*P[3][5]) + P[4][3]*(Dhat[0][3]*P[3][0] + Dhat[1][3]*P[3][1] + Dhat[2][3]*P[3][2] + Dhat[3][3]*P[3][3] + Dhat[4][3]*P[3][4] + Dhat[5][3]*P[3][5]) + P[4][4]*(Dhat[0][4]*P[3][0] + Dhat[1][4]*P[3][1] + Dhat[2][4]*P[3][2] + Dhat[3][4]*P[3][3] + Dhat[4][4]*P[3][4] + Dhat[5][4]*P[3][5]) + P[4][5]*(Dhat[0][5]*P[3][0] + Dhat[1][5]*P[3][1] + Dhat[2][5]*P[3][2] + Dhat[3][5]*P[3][3] + Dhat[4][5]*P[3][4] + Dhat[5][5]*P[3][5]); 
+  Diso[3][3] = P[3][0]*(Dhat[0][0]*P[3][0] + Dhat[1][0]*P[3][1] + Dhat[2][0]*P[3][2] + Dhat[3][0]*P[3][3] + Dhat[4][0]*P[3][4] + Dhat[5][0]*P[3][5]) + P[3][1]*(Dhat[0][1]*P[3][0] + Dhat[1][1]*P[3][1] + Dhat[2][1]*P[3][2] + Dhat[3][1]*P[3][3] + Dhat[4][1]*P[3][4] + Dhat[5][1]*P[3][5]) + P[3][2]*(Dhat[0][2]*P[3][0] + Dhat[1][2]*P[3][1] + Dhat[2][2]*P[3][2] + Dhat[3][2]*P[3][3] + Dhat[4][2]*P[3][4] + Dhat[5][2]*P[3][5]) + P[3][3]*(Dhat[0][3]*P[3][0] + Dhat[1][3]*P[3][1] + Dhat[2][3]*P[3][2] + Dhat[3][3]*P[3][3] + Dhat[4][3]*P[3][4] + Dhat[5][3]*P[3][5]) + P[3][4]*(Dhat[0][4]*P[3][0] + Dhat[1][4]*P[3][1] + Dhat[2][4]*P[3][2] + Dhat[3][4]*P[3][3] + Dhat[4][4]*P[3][4] + Dhat[5][4]*P[3][5]) + P[3][5]*(Dhat[0][5]*P[3][0] + Dhat[1][5]*P[3][1] + Dhat[2][5]*P[3][2] + Dhat[3][5]*P[3][3] + Dhat[4][5]*P[3][4] + Dhat[5][5]*P[3][5]);
+  Diso[3][4] = P[4][0]*(Dhat[0][0]*P[3][0] + Dhat[1][0]*P[3][1] + Dhat[2][0]*P[3][2] + Dhat[3][0]*P[3][3] + Dhat[4][0]*P[3][4] + Dhat[5][0]*P[3][5]) + P[4][1]*(Dhat[0][1]*P[3][0] + Dhat[1][1]*P[3][1] + Dhat[2][1]*P[3][2] + Dhat[3][1]*P[3][3] + Dhat[4][1]*P[3][4] + Dhat[5][1]*P[3][5]) + P[4][2]*(Dhat[0][2]*P[3][0] + Dhat[1][2]*P[3][1] + Dhat[2][2]*P[3][2] + Dhat[3][2]*P[3][3] + Dhat[4][2]*P[3][4] + Dhat[5][2]*P[3][5]) + P[4][3]*(Dhat[0][3]*P[3][0] + Dhat[1][3]*P[3][1] + Dhat[2][3]*P[3][2] + Dhat[3][3]*P[3][3] + Dhat[4][3]*P[3][4] + Dhat[5][3]*P[3][5]) + P[4][4]*(Dhat[0][4]*P[3][0] + Dhat[1][4]*P[3][1] + Dhat[2][4]*P[3][2] + Dhat[3][4]*P[3][3] + Dhat[4][4]*P[3][4] + Dhat[5][4]*P[3][5]) + P[4][5]*(Dhat[0][5]*P[3][0] + Dhat[1][5]*P[3][1] + Dhat[2][5]*P[3][2] + Dhat[3][5]*P[3][3] + Dhat[4][5]*P[3][4] + Dhat[5][5]*P[3][5]);
   Diso[3][5] = P[5][0]*(Dhat[0][0]*P[3][0] + Dhat[1][0]*P[3][1] + Dhat[2][0]*P[3][2] + Dhat[3][0]*P[3][3] + Dhat[4][0]*P[3][4] + Dhat[5][0]*P[3][5]) + P[5][1]*(Dhat[0][1]*P[3][0] + Dhat[1][1]*P[3][1] + Dhat[2][1]*P[3][2] + Dhat[3][1]*P[3][3] + Dhat[4][1]*P[3][4] + Dhat[5][1]*P[3][5]) + P[5][2]*(Dhat[0][2]*P[3][0] + Dhat[1][2]*P[3][1] + Dhat[2][2]*P[3][2] + Dhat[3][2]*P[3][3] + Dhat[4][2]*P[3][4] + Dhat[5][2]*P[3][5]) + P[5][3]*(Dhat[0][3]*P[3][0] + Dhat[1][3]*P[3][1] + Dhat[2][3]*P[3][2] + Dhat[3][3]*P[3][3] + Dhat[4][3]*P[3][4] + Dhat[5][3]*P[3][5]) + P[5][4]*(Dhat[0][4]*P[3][0] + Dhat[1][4]*P[3][1] + Dhat[2][4]*P[3][2] + Dhat[3][4]*P[3][3] + Dhat[4][4]*P[3][4] + Dhat[5][4]*P[3][5]) + P[5][5]*(Dhat[0][5]*P[3][0] + Dhat[1][5]*P[3][1] + Dhat[2][5]*P[3][2] + Dhat[3][5]*P[3][3] + Dhat[4][5]*P[3][4] + Dhat[5][5]*P[3][5]);
-  Diso[4][4] = P[4][0]*(Dhat[0][0]*P[4][0] + Dhat[1][0]*P[4][1] + Dhat[2][0]*P[4][2] + Dhat[3][0]*P[4][3] + Dhat[4][0]*P[4][4] + Dhat[5][0]*P[4][5]) + P[4][1]*(Dhat[0][1]*P[4][0] + Dhat[1][1]*P[4][1] + Dhat[2][1]*P[4][2] + Dhat[3][1]*P[4][3] + Dhat[4][1]*P[4][4] + Dhat[5][1]*P[4][5]) + P[4][2]*(Dhat[0][2]*P[4][0] + Dhat[1][2]*P[4][1] + Dhat[2][2]*P[4][2] + Dhat[3][2]*P[4][3] + Dhat[4][2]*P[4][4] + Dhat[5][2]*P[4][5]) + P[4][3]*(Dhat[0][3]*P[4][0] + Dhat[1][3]*P[4][1] + Dhat[2][3]*P[4][2] + Dhat[3][3]*P[4][3] + Dhat[4][3]*P[4][4] + Dhat[5][3]*P[4][5]) + P[4][4]*(Dhat[0][4]*P[4][0] + Dhat[1][4]*P[4][1] + Dhat[2][4]*P[4][2] + Dhat[3][4]*P[4][3] + Dhat[4][4]*P[4][4] + Dhat[5][4]*P[4][5]) + P[4][5]*(Dhat[0][5]*P[4][0] + Dhat[1][5]*P[4][1] + Dhat[2][5]*P[4][2] + Dhat[3][5]*P[4][3] + Dhat[4][5]*P[4][4] + Dhat[5][5]*P[4][5]); 
+  Diso[4][4] = P[4][0]*(Dhat[0][0]*P[4][0] + Dhat[1][0]*P[4][1] + Dhat[2][0]*P[4][2] + Dhat[3][0]*P[4][3] + Dhat[4][0]*P[4][4] + Dhat[5][0]*P[4][5]) + P[4][1]*(Dhat[0][1]*P[4][0] + Dhat[1][1]*P[4][1] + Dhat[2][1]*P[4][2] + Dhat[3][1]*P[4][3] + Dhat[4][1]*P[4][4] + Dhat[5][1]*P[4][5]) + P[4][2]*(Dhat[0][2]*P[4][0] + Dhat[1][2]*P[4][1] + Dhat[2][2]*P[4][2] + Dhat[3][2]*P[4][3] + Dhat[4][2]*P[4][4] + Dhat[5][2]*P[4][5]) + P[4][3]*(Dhat[0][3]*P[4][0] + Dhat[1][3]*P[4][1] + Dhat[2][3]*P[4][2] + Dhat[3][3]*P[4][3] + Dhat[4][3]*P[4][4] + Dhat[5][3]*P[4][5]) + P[4][4]*(Dhat[0][4]*P[4][0] + Dhat[1][4]*P[4][1] + Dhat[2][4]*P[4][2] + Dhat[3][4]*P[4][3] + Dhat[4][4]*P[4][4] + Dhat[5][4]*P[4][5]) + P[4][5]*(Dhat[0][5]*P[4][0] + Dhat[1][5]*P[4][1] + Dhat[2][5]*P[4][2] + Dhat[3][5]*P[4][3] + Dhat[4][5]*P[4][4] + Dhat[5][5]*P[4][5]);
   Diso[4][5] = P[5][0]*(Dhat[0][0]*P[4][0] + Dhat[1][0]*P[4][1] + Dhat[2][0]*P[4][2] + Dhat[3][0]*P[4][3] + Dhat[4][0]*P[4][4] + Dhat[5][0]*P[4][5]) + P[5][1]*(Dhat[0][1]*P[4][0] + Dhat[1][1]*P[4][1] + Dhat[2][1]*P[4][2] + Dhat[3][1]*P[4][3] + Dhat[4][1]*P[4][4] + Dhat[5][1]*P[4][5]) + P[5][2]*(Dhat[0][2]*P[4][0] + Dhat[1][2]*P[4][1] + Dhat[2][2]*P[4][2] + Dhat[3][2]*P[4][3] + Dhat[4][2]*P[4][4] + Dhat[5][2]*P[4][5]) + P[5][3]*(Dhat[0][3]*P[4][0] + Dhat[1][3]*P[4][1] + Dhat[2][3]*P[4][2] + Dhat[3][3]*P[4][3] + Dhat[4][3]*P[4][4] + Dhat[5][3]*P[4][5]) + P[5][4]*(Dhat[0][4]*P[4][0] + Dhat[1][4]*P[4][1] + Dhat[2][4]*P[4][2] + Dhat[3][4]*P[4][3] + Dhat[4][4]*P[4][4] + Dhat[5][4]*P[4][5]) + P[5][5]*(Dhat[0][5]*P[4][0] + Dhat[1][5]*P[4][1] + Dhat[2][5]*P[4][2] + Dhat[3][5]*P[4][3] + Dhat[4][5]*P[4][4] + Dhat[5][5]*P[4][5]);
   Diso[5][5] = P[5][0]*(Dhat[0][0]*P[5][0] + Dhat[1][0]*P[5][1] + Dhat[2][0]*P[5][2] + Dhat[3][0]*P[5][3] + Dhat[4][0]*P[5][4] + Dhat[5][0]*P[5][5]) + P[5][1]*(Dhat[0][1]*P[5][0] + Dhat[1][1]*P[5][1] + Dhat[2][1]*P[5][2] + Dhat[3][1]*P[5][3] + Dhat[4][1]*P[5][4] + Dhat[5][1]*P[5][5]) + P[5][2]*(Dhat[0][2]*P[5][0] + Dhat[1][2]*P[5][1] + Dhat[2][2]*P[5][2] + Dhat[3][2]*P[5][3] + Dhat[4][2]*P[5][4] + Dhat[5][2]*P[5][5]) + P[5][3]*(Dhat[0][3]*P[5][0] + Dhat[1][3]*P[5][1] + Dhat[2][3]*P[5][2] + Dhat[3][3]*P[5][3] + Dhat[4][3]*P[5][4] + Dhat[5][3]*P[5][5]) + P[5][4]*(Dhat[0][4]*P[5][0] + Dhat[1][4]*P[5][1] + Dhat[2][4]*P[5][2] + Dhat[3][4]*P[5][3] + Dhat[4][4]*P[5][4] + Dhat[5][4]*P[5][5]) + P[5][5]*(Dhat[0][5]*P[5][0] + Dhat[1][5]*P[5][1] + Dhat[2][5]*P[5][2] + Dhat[3][5]*P[5][3] + Dhat[4][5]*P[5][4] + Dhat[5][5]*P[5][5]);
-  
+
   Diso[0][0] += (2/3)*tempX*Pbar[0][0]-(2/3)*(Cinv[0][0]*Siso[0]+Siso[0]*Cinv[0][0]);
   Diso[0][1] += (2/3)*tempX*Pbar[0][1]-(2/3)*(Cinv[0][0]*Siso[1]+Siso[0]*Cinv[1][1]);
   Diso[0][2] += (2/3)*tempX*Pbar[0][2]-(2/3)*(Cinv[0][0]*Siso[2]+Siso[0]*Cinv[2][2]);
 #undef __FUNCT__
 #define __FUNCT__ "DeltaE"
 void DeltaE(PetscScalar Nx, PetscScalar Ny, PetscScalar Nz, PetscScalar (*F)[3], PetscScalar (*B)[3])
-{    
-  // Given F and basis values, returns B 
-  B[0][0] = F[0][0]*Nx; 
-  B[0][1] = F[1][0]*Nx; 
+{
+  // Given F and basis values, returns B
+  B[0][0] = F[0][0]*Nx;
+  B[0][1] = F[1][0]*Nx;
   B[0][2] = F[2][0]*Nx;
-  B[1][0] = F[0][1]*Ny; 
-  B[1][1] = F[1][1]*Ny; 
+  B[1][0] = F[0][1]*Ny;
+  B[1][1] = F[1][1]*Ny;
   B[1][2] = F[2][1]*Ny;
-  B[2][0] = F[0][2]*Nz; 
-  B[2][1] = F[1][2]*Nz; 
+  B[2][0] = F[0][2]*Nz;
+  B[2][1] = F[1][2]*Nz;
   B[2][2] = F[2][2]*Nz;
-  B[3][0] = F[0][0]*Ny+F[0][1]*Nx; 
-  B[3][1] = F[1][0]*Ny+F[1][1]*Nx; 
+  B[3][0] = F[0][0]*Ny+F[0][1]*Nx;
+  B[3][1] = F[1][0]*Ny+F[1][1]*Nx;
   B[3][2] = F[2][0]*Ny+F[2][1]*Nx;
-  B[4][0] = F[0][1]*Nz+F[0][2]*Ny; 
-  B[4][1] = F[1][1]*Nz+F[1][2]*Ny; 
+  B[4][0] = F[0][1]*Nz+F[0][2]*Ny;
+  B[4][1] = F[1][1]*Nz+F[1][2]*Ny;
   B[4][2] = F[2][1]*Nz+F[2][2]*Ny;
-  B[5][0] = F[0][0]*Nz+F[0][2]*Nx; 
-  B[5][1] = F[1][0]*Nz+F[1][2]*Nx; 
+  B[5][0] = F[0][0]*Nz+F[0][2]*Nx;
+  B[5][1] = F[1][0]*Nz+F[1][2]*Nx;
   B[5][2] = F[2][0]*Nz+F[2][2]*Nx;
 }
 
 #undef  __FUNCT__
 #define __FUNCT__ "Residual"
 PetscErrorCode Residual(IGAPoint pnt,const PetscScalar *U,PetscScalar *Re,void *ctx)
-{    
+{
   AppCtx *user = (AppCtx *)ctx;
 
   // call user model
     PetscReal Na_y  = N1[a][1];
     PetscReal Na_z  = N1[a][2];
     DeltaE(Na_x,Na_y,Na_z,F,B);
-    R[a][0] = B[0][0]*S[0][0]+B[1][0]*S[1][1]+B[2][0]*S[2][2]+B[3][0]*S[0][1]+B[4][0]*S[1][2]+B[5][0]*S[0][2]; 
+    R[a][0] = B[0][0]*S[0][0]+B[1][0]*S[1][1]+B[2][0]*S[2][2]+B[3][0]*S[0][1]+B[4][0]*S[1][2]+B[5][0]*S[0][2];
     R[a][1] = B[0][1]*S[0][0]+B[1][1]*S[1][1]+B[2][1]*S[2][2]+B[3][1]*S[0][1]+B[4][1]*S[1][2]+B[5][1]*S[0][2];
     R[a][2] = B[0][2]*S[0][0]+B[1][2]*S[1][1]+B[2][2]*S[2][2]+B[3][2]*S[0][1]+B[4][2]*S[1][2]+B[5][2]*S[0][2];
   }
 #undef  __FUNCT__
 #define __FUNCT__ "Jacobian"
 PetscErrorCode Jacobian(IGAPoint pnt,const PetscScalar *U,PetscScalar *Je,void *ctx)
-{    
+{
   AppCtx *user = (AppCtx *)ctx;
 
   // Call user model
 
     PetscReal Na_x  = N1[a][0];
     PetscReal Na_y  = N1[a][1];
-    PetscReal Na_z  = N1[a][2];  
+    PetscReal Na_z  = N1[a][2];
     DeltaE(Na_x,Na_y,Na_z,F,Ba);
 
     for (b=0; b<nen; b++) {
       PetscReal Nb_y  = N1[b][1];
       PetscReal Nb_z  = N1[b][2];
       DeltaE(Nb_x,Nb_y,Nb_z,F,Bb);
-      
-      Chi[0][0]=Bb[0][0]*(Ba[0][0]*D[0][0] + Ba[1][0]*D[1][0] + Ba[2][0]*D[2][0] + Ba[3][0]*D[3][0] + Ba[4][0]*D[4][0] + Ba[5][0]*D[5][0]) + 
-        Bb[1][0]*(Ba[0][0]*D[0][1] + Ba[1][0]*D[1][1] + Ba[2][0]*D[2][1] + Ba[3][0]*D[3][1] + Ba[4][0]*D[4][1] + Ba[5][0]*D[5][1]) + 
-        Bb[2][0]*(Ba[0][0]*D[0][2] + Ba[1][0]*D[1][2] + Ba[2][0]*D[2][2] + Ba[3][0]*D[3][2] + Ba[4][0]*D[4][2] + Ba[5][0]*D[5][2]) + 
-        Bb[3][0]*(Ba[0][0]*D[0][3] + Ba[1][0]*D[1][3] + Ba[2][0]*D[2][3] + Ba[3][0]*D[3][3] + Ba[4][0]*D[4][3] + Ba[5][0]*D[5][3]) + 
-        Bb[4][0]*(Ba[0][0]*D[0][4] + Ba[1][0]*D[1][4] + Ba[2][0]*D[2][4] + Ba[3][0]*D[3][4] + Ba[4][0]*D[4][4] + Ba[5][0]*D[5][4]) + 
+
+      Chi[0][0]=Bb[0][0]*(Ba[0][0]*D[0][0] + Ba[1][0]*D[1][0] + Ba[2][0]*D[2][0] + Ba[3][0]*D[3][0] + Ba[4][0]*D[4][0] + Ba[5][0]*D[5][0]) +
+        Bb[1][0]*(Ba[0][0]*D[0][1] + Ba[1][0]*D[1][1] + Ba[2][0]*D[2][1] + Ba[3][0]*D[3][1] + Ba[4][0]*D[4][1] + Ba[5][0]*D[5][1]) +
+        Bb[2][0]*(Ba[0][0]*D[0][2] + Ba[1][0]*D[1][2] + Ba[2][0]*D[2][2] + Ba[3][0]*D[3][2] + Ba[4][0]*D[4][2] + Ba[5][0]*D[5][2]) +
+        Bb[3][0]*(Ba[0][0]*D[0][3] + Ba[1][0]*D[1][3] + Ba[2][0]*D[2][3] + Ba[3][0]*D[3][3] + Ba[4][0]*D[4][3] + Ba[5][0]*D[5][3]) +
+        Bb[4][0]*(Ba[0][0]*D[0][4] + Ba[1][0]*D[1][4] + Ba[2][0]*D[2][4] + Ba[3][0]*D[3][4] + Ba[4][0]*D[4][4] + Ba[5][0]*D[5][4]) +
         Bb[5][0]*(Ba[0][0]*D[0][5] + Ba[1][0]*D[1][5] + Ba[2][0]*D[2][5] + Ba[3][0]*D[3][5] + Ba[4][0]*D[4][5] + Ba[5][0]*D[5][5]);
-        
-      Chi[0][1]=Bb[0][1]*(Ba[0][0]*D[0][0] + Ba[1][0]*D[1][0] + Ba[2][0]*D[2][0] + Ba[3][0]*D[3][0] + Ba[4][0]*D[4][0] + Ba[5][0]*D[5][0]) + 
-        Bb[1][1]*(Ba[0][0]*D[0][1] + Ba[1][0]*D[1][1] + Ba[2][0]*D[2][1] + Ba[3][0]*D[3][1] + Ba[4][0]*D[4][1] + Ba[5][0]*D[5][1]) + 
-        Bb[2][1]*(Ba[0][0]*D[0][2] + Ba[1][0]*D[1][2] + Ba[2][0]*D[2][2] + Ba[3][0]*D[3][2] + Ba[4][0]*D[4][2] + Ba[5][0]*D[5][2]) + 
-        Bb[3][1]*(Ba[0][0]*D[0][3] + Ba[1][0]*D[1][3] + Ba[2][0]*D[2][3] + Ba[3][0]*D[3][3] + Ba[4][0]*D[4][3] + Ba[5][0]*D[5][3]) + 
-        Bb[4][1]*(Ba[0][0]*D[0][4] + Ba[1][0]*D[1][4] + Ba[2][0]*D[2][4] + Ba[3][0]*D[3][4] + Ba[4][0]*D[4][4] + Ba[5][0]*D[5][4]) + 
+
+      Chi[0][1]=Bb[0][1]*(Ba[0][0]*D[0][0] + Ba[1][0]*D[1][0] + Ba[2][0]*D[2][0] + Ba[3][0]*D[3][0] + Ba[4][0]*D[4][0] + Ba[5][0]*D[5][0]) +
+        Bb[1][1]*(Ba[0][0]*D[0][1] + Ba[1][0]*D[1][1] + Ba[2][0]*D[2][1] + Ba[3][0]*D[3][1] + Ba[4][0]*D[4][1] + Ba[5][0]*D[5][1]) +
+        Bb[2][1]*(Ba[0][0]*D[0][2] + Ba[1][0]*D[1][2] + Ba[2][0]*D[2][2] + Ba[3][0]*D[3][2] + Ba[4][0]*D[4][2] + Ba[5][0]*D[5][2]) +
+        Bb[3][1]*(Ba[0][0]*D[0][3] + Ba[1][0]*D[1][3] + Ba[2][0]*D[2][3] + Ba[3][0]*D[3][3] + Ba[4][0]*D[4][3] + Ba[5][0]*D[5][3]) +
+        Bb[4][1]*(Ba[0][0]*D[0][4] + Ba[1][0]*D[1][4] + Ba[2][0]*D[2][4] + Ba[3][0]*D[3][4] + Ba[4][0]*D[4][4] + Ba[5][0]*D[5][4]) +
         Bb[5][1]*(Ba[0][0]*D[0][5] + Ba[1][0]*D[1][5] + Ba[2][0]*D[2][5] + Ba[3][0]*D[3][5] + Ba[4][0]*D[4][5] + Ba[5][0]*D[5][5]);
-        
-      Chi[0][2]=Bb[0][2]*(Ba[0][0]*D[0][0] + Ba[1][0]*D[1][0] + Ba[2][0]*D[2][0] + Ba[3][0]*D[3][0] + Ba[4][0]*D[4][0] + Ba[5][0]*D[5][0]) + 
-        Bb[1][2]*(Ba[0][0]*D[0][1] + Ba[1][0]*D[1][1] + Ba[2][0]*D[2][1] + Ba[3][0]*D[3][1] + Ba[4][0]*D[4][1] + Ba[5][0]*D[5][1]) + 
-        Bb[2][2]*(Ba[0][0]*D[0][2] + Ba[1][0]*D[1][2] + Ba[2][0]*D[2][2] + Ba[3][0]*D[3][2] + Ba[4][0]*D[4][2] + Ba[5][0]*D[5][2]) + 
-        Bb[3][2]*(Ba[0][0]*D[0][3] + Ba[1][0]*D[1][3] + Ba[2][0]*D[2][3] + Ba[3][0]*D[3][3] + Ba[4][0]*D[4][3] + Ba[5][0]*D[5][3]) + 
-        Bb[4][2]*(Ba[0][0]*D[0][4] + Ba[1][0]*D[1][4] + Ba[2][0]*D[2][4] + Ba[3][0]*D[3][4] + Ba[4][0]*D[4][4] + Ba[5][0]*D[5][4]) + 
+
+      Chi[0][2]=Bb[0][2]*(Ba[0][0]*D[0][0] + Ba[1][0]*D[1][0] + Ba[2][0]*D[2][0] + Ba[3][0]*D[3][0] + Ba[4][0]*D[4][0] + Ba[5][0]*D[5][0]) +
+        Bb[1][2]*(Ba[0][0]*D[0][1] + Ba[1][0]*D[1][1] + Ba[2][0]*D[2][1] + Ba[3][0]*D[3][1] + Ba[4][0]*D[4][1] + Ba[5][0]*D[5][1]) +
+        Bb[2][2]*(Ba[0][0]*D[0][2] + Ba[1][0]*D[1][2] + Ba[2][0]*D[2][2] + Ba[3][0]*D[3][2] + Ba[4][0]*D[4][2] + Ba[5][0]*D[5][2]) +
+        Bb[3][2]*(Ba[0][0]*D[0][3] + Ba[1][0]*D[1][3] + Ba[2][0]*D[2][3] + Ba[3][0]*D[3][3] + Ba[4][0]*D[4][3] + Ba[5][0]*D[5][3]) +
+        Bb[4][2]*(Ba[0][0]*D[0][4] + Ba[1][0]*D[1][4] + Ba[2][0]*D[2][4] + Ba[3][0]*D[3][4] + Ba[4][0]*D[4][4] + Ba[5][0]*D[5][4]) +
         Bb[5][2]*(Ba[0][0]*D[0][5] + Ba[1][0]*D[1][5] + Ba[2][0]*D[2][5] + Ba[3][0]*D[3][5] + Ba[4][0]*D[4][5] + Ba[5][0]*D[5][5]);
-        
-      Chi[1][0]=Bb[0][0]*(Ba[0][1]*D[0][0] + Ba[1][1]*D[1][0] + Ba[2][1]*D[2][0] + Ba[3][1]*D[3][0] + Ba[4][1]*D[4][0] + Ba[5][1]*D[5][0]) + 
-        Bb[1][0]*(Ba[0][1]*D[0][1] + Ba[1][1]*D[1][1] + Ba[2][1]*D[2][1] + Ba[3][1]*D[3][1] + Ba[4][1]*D[4][1] + Ba[5][1]*D[5][1]) + 
-        Bb[2][0]*(Ba[0][1]*D[0][2] + Ba[1][1]*D[1][2] + Ba[2][1]*D[2][2] + Ba[3][1]*D[3][2] + Ba[4][1]*D[4][2] + Ba[5][1]*D[5][2]) + 
-        Bb[3][0]*(Ba[0][1]*D[0][3] + Ba[1][1]*D[1][3] + Ba[2][1]*D[2][3] + Ba[3][1]*D[3][3] + Ba[4][1]*D[4][3] + Ba[5][1]*D[5][3]) + 
-        Bb[4][0]*(Ba[0][1]*D[0][4] + Ba[1][1]*D[1][4] + Ba[2][1]*D[2][4] + Ba[3][1]*D[3][4] + Ba[4][1]*D[4][4] + Ba[5][1]*D[5][4]) + 
+
+      Chi[1][0]=Bb[0][0]*(Ba[0][1]*D[0][0] + Ba[1][1]*D[1][0] + Ba[2][1]*D[2][0] + Ba[3][1]*D[3][0] + Ba[4][1]*D[4][0] + Ba[5][1]*D[5][0]) +
+        Bb[1][0]*(Ba[0][1]*D[0][1] + Ba[1][1]*D[1][1] + Ba[2][1]*D[2][1] + Ba[3][1]*D[3][1] + Ba[4][1]*D[4][1] + Ba[5][1]*D[5][1]) +
+        Bb[2][0]*(Ba[0][1]*D[0][2] + Ba[1][1]*D[1][2] + Ba[2][1]*D[2][2] + Ba[3][1]*D[3][2] + Ba[4][1]*D[4][2] + Ba[5][1]*D[5][2]) +
+        Bb[3][0]*(Ba[0][1]*D[0][3] + Ba[1][1]*D[1][3] + Ba[2][1]*D[2][3] + Ba[3][1]*D[3][3] + Ba[4][1]*D[4][3] + Ba[5][1]*D[5][3]) +
+        Bb[4][0]*(Ba[0][1]*D[0][4] + Ba[1][1]*D[1][4] + Ba[2][1]*D[2][4] + Ba[3][1]*D[3][4] + Ba[4][1]*D[4][4] + Ba[5][1]*D[5][4]) +
         Bb[5][0]*(Ba[0][1]*D[0][5] + Ba[1][1]*D[1][5] + Ba[2][1]*D[2][5] + Ba[3][1]*D[3][5] + Ba[4][1]*D[4][5] + Ba[5][1]*D[5][5]);
-        
-      Chi[1][1]=Bb[0][1]*(Ba[0][1]*D[0][0] + Ba[1][1]*D[1][0] + Ba[2][1]*D[2][0] + Ba[3][1]*D[3][0] + Ba[4][1]*D[4][0] + Ba[5][1]*D[5][0]) + 
-        Bb[1][1]*(Ba[0][1]*D[0][1] + Ba[1][1]*D[1][1] + Ba[2][1]*D[2][1] + Ba[3][1]*D[3][1] + Ba[4][1]*D[4][1] + Ba[5][1]*D[5][1]) + 
-        Bb[2][1]*(Ba[0][1]*D[0][2] + Ba[1][1]*D[1][2] + Ba[2][1]*D[2][2] + Ba[3][1]*D[3][2] + Ba[4][1]*D[4][2] + Ba[5][1]*D[5][2]) + 
-        Bb[3][1]*(Ba[0][1]*D[0][3] + Ba[1][1]*D[1][3] + Ba[2][1]*D[2][3] + Ba[3][1]*D[3][3] + Ba[4][1]*D[4][3] + Ba[5][1]*D[5][3]) + 
-        Bb[4][1]*(Ba[0][1]*D[0][4] + Ba[1][1]*D[1][4] + Ba[2][1]*D[2][4] + Ba[3][1]*D[3][4] + Ba[4][1]*D[4][4] + Ba[5][1]*D[5][4]) + 
+
+      Chi[1][1]=Bb[0][1]*(Ba[0][1]*D[0][0] + Ba[1][1]*D[1][0] + Ba[2][1]*D[2][0] + Ba[3][1]*D[3][0] + Ba[4][1]*D[4][0] + Ba[5][1]*D[5][0]) +
+        Bb[1][1]*(Ba[0][1]*D[0][1] + Ba[1][1]*D[1][1] + Ba[2][1]*D[2][1] + Ba[3][1]*D[3][1] + Ba[4][1]*D[4][1] + Ba[5][1]*D[5][1]) +
+        Bb[2][1]*(Ba[0][1]*D[0][2] + Ba[1][1]*D[1][2] + Ba[2][1]*D[2][2] + Ba[3][1]*D[3][2] + Ba[4][1]*D[4][2] + Ba[5][1]*D[5][2]) +
+        Bb[3][1]*(Ba[0][1]*D[0][3] + Ba[1][1]*D[1][3] + Ba[2][1]*D[2][3] + Ba[3][1]*D[3][3] + Ba[4][1]*D[4][3] + Ba[5][1]*D[5][3]) +
+        Bb[4][1]*(Ba[0][1]*D[0][4] + Ba[1][1]*D[1][4] + Ba[2][1]*D[2][4] + Ba[3][1]*D[3][4] + Ba[4][1]*D[4][4] + Ba[5][1]*D[5][4]) +
         Bb[5][1]*(Ba[0][1]*D[0][5] + Ba[1][1]*D[1][5] + Ba[2][1]*D[2][5] + Ba[3][1]*D[3][5] + Ba[4][1]*D[4][5] + Ba[5][1]*D[5][5]);
-        
-      Chi[1][2]=Bb[0][2]*(Ba[0][1]*D[0][0] + Ba[1][1]*D[1][0] + Ba[2][1]*D[2][0] + Ba[3][1]*D[3][0] + Ba[4][1]*D[4][0] + Ba[5][1]*D[5][0]) + 
-        Bb[1][2]*(Ba[0][1]*D[0][1] + Ba[1][1]*D[1][1] + Ba[2][1]*D[2][1] + Ba[3][1]*D[3][1] + Ba[4][1]*D[4][1] + Ba[5][1]*D[5][1]) + 
-        Bb[2][2]*(Ba[0][1]*D[0][2] + Ba[1][1]*D[1][2] + Ba[2][1]*D[2][2] + Ba[3][1]*D[3][2] + Ba[4][1]*D[4][2] + Ba[5][1]*D[5][2]) + 
-        Bb[3][2]*(Ba[0][1]*D[0][3] + Ba[1][1]*D[1][3] + Ba[2][1]*D[2][3] + Ba[3][1]*D[3][3] + Ba[4][1]*D[4][3] + Ba[5][1]*D[5][3]) + 
-        Bb[4][2]*(Ba[0][1]*D[0][4] + Ba[1][1]*D[1][4] + Ba[2][1]*D[2][4] + Ba[3][1]*D[3][4] + Ba[4][1]*D[4][4] + Ba[5][1]*D[5][4]) + 
+
+      Chi[1][2]=Bb[0][2]*(Ba[0][1]*D[0][0] + Ba[1][1]*D[1][0] + Ba[2][1]*D[2][0] + Ba[3][1]*D[3][0] + Ba[4][1]*D[4][0] + Ba[5][1]*D[5][0]) +
+        Bb[1][2]*(Ba[0][1]*D[0][1] + Ba[1][1]*D[1][1] + Ba[2][1]*D[2][1] + Ba[3][1]*D[3][1] + Ba[4][1]*D[4][1] + Ba[5][1]*D[5][1]) +
+        Bb[2][2]*(Ba[0][1]*D[0][2] + Ba[1][1]*D[1][2] + Ba[2][1]*D[2][2] + Ba[3][1]*D[3][2] + Ba[4][1]*D[4][2] + Ba[5][1]*D[5][2]) +
+        Bb[3][2]*(Ba[0][1]*D[0][3] + Ba[1][1]*D[1][3] + Ba[2][1]*D[2][3] + Ba[3][1]*D[3][3] + Ba[4][1]*D[4][3] + Ba[5][1]*D[5][3]) +
+        Bb[4][2]*(Ba[0][1]*D[0][4] + Ba[1][1]*D[1][4] + Ba[2][1]*D[2][4] + Ba[3][1]*D[3][4] + Ba[4][1]*D[4][4] + Ba[5][1]*D[5][4]) +
         Bb[5][2]*(Ba[0][1]*D[0][5] + Ba[1][1]*D[1][5] + Ba[2][1]*D[2][5] + Ba[3][1]*D[3][5] + Ba[4][1]*D[4][5] + Ba[5][1]*D[5][5]);
-        
-      Chi[2][0]=Bb[0][0]*(Ba[0][2]*D[0][0] + Ba[1][2]*D[1][0] + Ba[2][2]*D[2][0] + Ba[3][2]*D[3][0] + Ba[4][2]*D[4][0] + Ba[5][2]*D[5][0]) + 
-        Bb[1][0]*(Ba[0][2]*D[0][1] + Ba[1][2]*D[1][1] + Ba[2][2]*D[2][1] + Ba[3][2]*D[3][1] + Ba[4][2]*D[4][1] + Ba[5][2]*D[5][1]) + 
-        Bb[2][0]*(Ba[0][2]*D[0][2] + Ba[1][2]*D[1][2] + Ba[2][2]*D[2][2] + Ba[3][2]*D[3][2] + Ba[4][2]*D[4][2] + Ba[5][2]*D[5][2]) + 
-        Bb[3][0]*(Ba[0][2]*D[0][3] + Ba[1][2]*D[1][3] + Ba[2][2]*D[2][3] + Ba[3][2]*D[3][3] + Ba[4][2]*D[4][3] + Ba[5][2]*D[5][3]) + 
-        Bb[4][0]*(Ba[0][2]*D[0][4] + Ba[1][2]*D[1][4] + Ba[2][2]*D[2][4] + Ba[3][2]*D[3][4] + Ba[4][2]*D[4][4] + Ba[5][2]*D[5][4]) + 
-        Bb[5][0]*(Ba[0][2]*D[0][5] + Ba[1][2]*D[1][5] + Ba[2][2]*D[2][5] + Ba[3][2]*D[3][5] + Ba[4][2]*D[4][5] + Ba[5][2]*D[5][5]); 
-        
-      Chi[2][1]=Bb[0][1]*(Ba[0][2]*D[0][0] + Ba[1][2]*D[1][0] + Ba[2][2]*D[2][0] + Ba[3][2]*D[3][0] + Ba[4][2]*D[4][0] + Ba[5][2]*D[5][0]) + 
-        Bb[1][1]*(Ba[0][2]*D[0][1] + Ba[1][2]*D[1][1] + Ba[2][2]*D[2][1] + Ba[3][2]*D[3][1] + Ba[4][2]*D[4][1] + Ba[5][2]*D[5][1]) + 
-        Bb[2][1]*(Ba[0][2]*D[0][2] + Ba[1][2]*D[1][2] + Ba[2][2]*D[2][2] + Ba[3][2]*D[3][2] + Ba[4][2]*D[4][2] + Ba[5][2]*D[5][2]) + 
-        Bb[3][1]*(Ba[0][2]*D[0][3] + Ba[1][2]*D[1][3] + Ba[2][2]*D[2][3] + Ba[3][2]*D[3][3] + Ba[4][2]*D[4][3] + Ba[5][2]*D[5][3]) + 
-        Bb[4][1]*(Ba[0][2]*D[0][4] + Ba[1][2]*D[1][4] + Ba[2][2]*D[2][4] + Ba[3][2]*D[3][4] + Ba[4][2]*D[4][4] + Ba[5][2]*D[5][4]) + 
-        Bb[5][1]*(Ba[0][2]*D[0][5] + Ba[1][2]*D[1][5] + Ba[2][2]*D[2][5] + Ba[3][2]*D[3][5] + Ba[4][2]*D[4][5] + Ba[5][2]*D[5][5]); 
-        
-      Chi[2][2]=Bb[0][2]*(Ba[0][2]*D[0][0] + Ba[1][2]*D[1][0] + Ba[2][2]*D[2][0] + Ba[3][2]*D[3][0] + Ba[4][2]*D[4][0] + Ba[5][2]*D[5][0]) + 
-        Bb[1][2]*(Ba[0][2]*D[0][1] + Ba[1][2]*D[1][1] + Ba[2][2]*D[2][1] + Ba[3][2]*D[3][1] + Ba[4][2]*D[4][1] + Ba[5][2]*D[5][1]) + 
-        Bb[2][2]*(Ba[0][2]*D[0][2] + Ba[1][2]*D[1][2] + Ba[2][2]*D[2][2] + Ba[3][2]*D[3][2] + Ba[4][2]*D[4][2] + Ba[5][2]*D[5][2]) + 
-        Bb[3][2]*(Ba[0][2]*D[0][3] + Ba[1][2]*D[1][3] + Ba[2][2]*D[2][3] + Ba[3][2]*D[3][3] + Ba[4][2]*D[4][3] + Ba[5][2]*D[5][3]) + 
-        Bb[4][2]*(Ba[0][2]*D[0][4] + Ba[1][2]*D[1][4] + Ba[2][2]*D[2][4] + Ba[3][2]*D[3][4] + Ba[4][2]*D[4][4] + Ba[5][2]*D[5][4]) + 
+
+      Chi[2][0]=Bb[0][0]*(Ba[0][2]*D[0][0] + Ba[1][2]*D[1][0] + Ba[2][2]*D[2][0] + Ba[3][2]*D[3][0] + Ba[4][2]*D[4][0] + Ba[5][2]*D[5][0]) +
+        Bb[1][0]*(Ba[0][2]*D[0][1] + Ba[1][2]*D[1][1] + Ba[2][2]*D[2][1] + Ba[3][2]*D[3][1] + Ba[4][2]*D[4][1] + Ba[5][2]*D[5][1]) +
+        Bb[2][0]*(Ba[0][2]*D[0][2] + Ba[1][2]*D[1][2] + Ba[2][2]*D[2][2] + Ba[3][2]*D[3][2] + Ba[4][2]*D[4][2] + Ba[5][2]*D[5][2]) +
+        Bb[3][0]*(Ba[0][2]*D[0][3] + Ba[1][2]*D[1][3] + Ba[2][2]*D[2][3] + Ba[3][2]*D[3][3] + Ba[4][2]*D[4][3] + Ba[5][2]*D[5][3]) +
+        Bb[4][0]*(Ba[0][2]*D[0][4] + Ba[1][2]*D[1][4] + Ba[2][2]*D[2][4] + Ba[3][2]*D[3][4] + Ba[4][2]*D[4][4] + Ba[5][2]*D[5][4]) +
+        Bb[5][0]*(Ba[0][2]*D[0][5] + Ba[1][2]*D[1][5] + Ba[2][2]*D[2][5] + Ba[3][2]*D[3][5] + Ba[4][2]*D[4][5] + Ba[5][2]*D[5][5]);
+
+      Chi[2][1]=Bb[0][1]*(Ba[0][2]*D[0][0] + Ba[1][2]*D[1][0] + Ba[2][2]*D[2][0] + Ba[3][2]*D[3][0] + Ba[4][2]*D[4][0] + Ba[5][2]*D[5][0]) +
+        Bb[1][1]*(Ba[0][2]*D[0][1] + Ba[1][2]*D[1][1] + Ba[2][2]*D[2][1] + Ba[3][2]*D[3][1] + Ba[4][2]*D[4][1] + Ba[5][2]*D[5][1]) +
+        Bb[2][1]*(Ba[0][2]*D[0][2] + Ba[1][2]*D[1][2] + Ba[2][2]*D[2][2] + Ba[3][2]*D[3][2] + Ba[4][2]*D[4][2] + Ba[5][2]*D[5][2]) +
+        Bb[3][1]*(Ba[0][2]*D[0][3] + Ba[1][2]*D[1][3] + Ba[2][2]*D[2][3] + Ba[3][2]*D[3][3] + Ba[4][2]*D[4][3] + Ba[5][2]*D[5][3]) +
+        Bb[4][1]*(Ba[0][2]*D[0][4] + Ba[1][2]*D[1][4] + Ba[2][2]*D[2][4] + Ba[3][2]*D[3][4] + Ba[4][2]*D[4][4] + Ba[5][2]*D[5][4]) +
+        Bb[5][1]*(Ba[0][2]*D[0][5] + Ba[1][2]*D[1][5] + Ba[2][2]*D[2][5] + Ba[3][2]*D[3][5] + Ba[4][2]*D[4][5] + Ba[5][2]*D[5][5]);
+
+      Chi[2][2]=Bb[0][2]*(Ba[0][2]*D[0][0] + Ba[1][2]*D[1][0] + Ba[2][2]*D[2][0] + Ba[3][2]*D[3][0] + Ba[4][2]*D[4][0] + Ba[5][2]*D[5][0]) +
+        Bb[1][2]*(Ba[0][2]*D[0][1] + Ba[1][2]*D[1][1] + Ba[2][2]*D[2][1] + Ba[3][2]*D[3][1] + Ba[4][2]*D[4][1] + Ba[5][2]*D[5][1]) +
+        Bb[2][2]*(Ba[0][2]*D[0][2] + Ba[1][2]*D[1][2] + Ba[2][2]*D[2][2] + Ba[3][2]*D[3][2] + Ba[4][2]*D[4][2] + Ba[5][2]*D[5][2]) +
+        Bb[3][2]*(Ba[0][2]*D[0][3] + Ba[1][2]*D[1][3] + Ba[2][2]*D[2][3] + Ba[3][2]*D[3][3] + Ba[4][2]*D[4][3] + Ba[5][2]*D[5][3]) +
+        Bb[4][2]*(Ba[0][2]*D[0][4] + Ba[1][2]*D[1][4] + Ba[2][2]*D[2][4] + Ba[3][2]*D[3][4] + Ba[4][2]*D[4][4] + Ba[5][2]*D[5][4]) +
         Bb[5][2]*(Ba[0][2]*D[0][5] + Ba[1][2]*D[1][5] + Ba[2][2]*D[2][5] + Ba[3][2]*D[3][5] + Ba[4][2]*D[4][5] + Ba[5][2]*D[5][5]);
 
       G=Na_x*(S[0][0]*Nb_x + S[0][1]*Nb_y + S[0][2]*Nb_z) +
         Na_y*(S[1][0]*Nb_x + S[1][1]*Nb_y + S[1][2]*Nb_z) +
         Na_z*(S[2][0]*Nb_x + S[2][1]*Nb_y + S[2][2]*Nb_z);
-      
+
       K[a][0][b][0] = G+Chi[0][0];
       K[a][1][b][0] =   Chi[1][0];
       K[a][2][b][0] =   Chi[2][0];
-      
+
       K[a][0][b][1] =   Chi[0][1];
       K[a][1][b][1] = G+Chi[1][1];
       K[a][2][b][1] =   Chi[2][1];
-      
+
       K[a][0][b][2] =   Chi[0][2];
       K[a][1][b][2] =   Chi[1][2];
       K[a][2][b][2] = G+Chi[2][2];
     }
   }
-  
+
   return 0;
 }
 
 
 #undef __FUNCT__
 #define __FUNCT__ "main"
-int main(int argc, char *argv[]) 
+int main(int argc, char *argv[])
 {
   // Initialization of PETSc
   PetscErrorCode ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
-  
+
   // Application specific data (defaults to Aluminum)
   AppCtx user;
-  PetscScalar E  = 70.0e9; 
+  PetscScalar E  = 70.0e9;
   PetscScalar nu = 0.35;
   user.model     = GeneralModel;
   PetscBool NeoHook,StVenant,MooneyR1,MooneyR2;
   NeoHook = StVenant = MooneyR1 = MooneyR2 = PETSC_FALSE;
   PetscInt nsteps = 1;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","HyperElasticity Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-neohook","Use the NeoHookean constitutive model",__FILE__,NeoHook,&NeoHook,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-mooneyr1","Use the MooneyRivlin1 constitutive model",__FILE__,MooneyR1,&MooneyR1,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-mooneyr2","Use the MooneyRivlin2 constitutive model",__FILE__,MooneyR2,&MooneyR2,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-stvenant","Use the StVenant constitutive model",__FILE__,StVenant,&StVenant,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-nsteps","Number of load steps to take",__FILE__,nsteps,&nsteps,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-neohook","Use the NeoHookean constitutive model",__FILE__,NeoHook,&NeoHook,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-mooneyr1","Use the MooneyRivlin1 constitutive model",__FILE__,MooneyR1,&MooneyR1,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-mooneyr2","Use the MooneyRivlin2 constitutive model",__FILE__,MooneyR2,&MooneyR2,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-stvenant","Use the StVenant constitutive model",__FILE__,StVenant,&StVenant,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-nsteps","Number of load steps to take",__FILE__,nsteps,&nsteps,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
-  
+
   user.lambda = E*nu/(1.+nu)/(1.-2.*nu);
   user.mu     = 0.5*E/(1.+nu);
   user.c1     = user.mu*0.5;
   // Load stepping
   PetscInt step;
   for(step=0;step<nsteps;step++){
-    
+
     PetscPrintf(PETSC_COMM_WORLD,"%d Load Step\n",step);
-    
+
     // Solve step
     ierr = VecZeroEntries(U);CHKERRQ(ierr);
-    ierr = SNESSolve(snes,PETSC_NULL,U);CHKERRQ(ierr);
+    ierr = SNESSolve(snes,NULL,U);CHKERRQ(ierr);
 
     // Store total displacement
     ierr = VecAXPY(Utotal,1.0,U);CHKERRQ(ierr);
   return (x-1)*(x-a)*(x+a)*(x+1);
 }
 
+PetscScalar Hill(PetscReal x)
+{
+  return (x-1)*(x-1)*(x+1)*(x+1);
+}
+
+PetscScalar Sine(PetscReal x)
+{
+  return sin(M_PI*x);
+}
+
 
 typedef struct {
   PetscScalar (*Function)(PetscReal x);
 {
   AppCtx *user = (AppCtx *)ctx;
   PetscReal x = p->point[0];
+  PetscScalar f = user->Function(x);
+
   PetscReal *N = p->shape[0];
-  
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na = N[a];
       PetscReal Nb = N[b];
       K[a*nen+b] = Na * Nb;
     }
-    F[a] = Na * user->Function(x);
+    F[a] = Na * f;
   }
   return 0;
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal x = p->point[0];
+  PetscScalar f = user->Function(x);
+
+  PetscScalar u;
+  IGAPointFormValue(p,U,&u);
+
+  PetscReal e = PetscAbsScalar(u - f);
+  S[0] = e*e;
+
+  return 0;
+}
+
 #undef __FUNCT__
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
   AppCtx user;
 
-  PetscInt choice=2;
-  const char *choicelist[] = {"line", "parabola", "poly3", "poly4", 0};
-  PetscBool t=PETSC_FALSE;
-  PetscInt N=16, p=2, C=PETSC_DECIDE;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Projection1D Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-periodic", "periodic",__FILE__,t,&t,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-N", "number of elements (along one dimension)",__FILE__,N,&N,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-p", "polynomial order",__FILE__,p,&p,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-C", "global continuity order",__FILE__,C,&C,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsEList("-function","1D function",__FILE__,choicelist,4,choicelist[choice],&choice,PETSC_NULL);CHKERRQ(ierr);
+  PetscBool w=PETSC_FALSE;
+  PetscInt  N=16;
+  PetscInt  p=2;
+  PetscInt  C=PETSC_DECIDE;
+  PetscInt  choice=2;
+  const char *choicelist[] = {"line", "parabola", "poly3", "poly4", "hill", "sine", 0};
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","L2 Projection 1D Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBool ("-periodic", "periodicity",__FILE__,w,&w,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt  ("-N", "number of elements",__FILE__,N,&N,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt  ("-p", "polynomial order",  __FILE__,p,&p,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt  ("-C", "continuity order",  __FILE__,C,&C,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEList("-function","1D function", __FILE__,choicelist,6,choicelist[choice],&choice,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
   if (C == PETSC_DECIDE) C = p-1;
   switch (choice) {
   case 1: user.Function = Parabola; break;
   case 2: user.Function = Poly3;    break;
   case 3: user.Function = Poly4;    break;
+  case 4: user.Function = Hill;     break;
+  case 5: user.Function = Sine;     break;
   }
 
   IGA iga;
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
   IGAAxis axis;
   ierr = IGAGetAxis(iga,0,&axis);CHKERRQ(ierr);
-  ierr = IGAAxisSetPeriodic(axis,t);CHKERRQ(ierr);
+  ierr = IGAAxisSetPeriodic(axis,w);CHKERRQ(ierr);
   ierr = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
   ierr = IGAAxisInitUniform(axis,N,-1.0,1.0,C);CHKERRQ(ierr);
 
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
   ierr = IGASetUserSystem(iga,System,&user);CHKERRQ(ierr);
   ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
-  
+
   KSP ksp;
   ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
   ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
-  ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
-  
+  PetscScalar error = 0;
+  ierr = IGAFormScalar(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);
+  if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);}
+
+  PetscBool draw = PETSC_FALSE;
+  ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
+  if (draw) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);
     /**/ - 1.0/3 * exp(-pow(X+1,2) - pow(Y,2));
 }
 
+PetscScalar Hill(PetscReal x, PetscReal y)
+{
+  return   (x-1)*(x-1)*(x+1)*(x+1)
+    /**/ * (y-1)*(y-1)*(y+1)*(y+1);
+}
+
+PetscScalar Sine(PetscReal x, PetscReal y)
+{
+  return   sin(M_PI*x)
+    /**/ * sin(M_PI*y);
+}
+
 typedef struct {
   PetscScalar (*Function)(PetscReal x, PetscReal y);
 } AppCtx;
   AppCtx *user = (AppCtx *)ctx;
   PetscReal x = p->point[0];
   PetscReal y = p->point[1];
+  PetscScalar f = user->Function(x,y);
+
   PetscReal *N = p->shape[0];
-
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na = N[a];
       PetscReal Nb = N[b];
       K[a*nen+b] = Na * Nb;
     }
-    F[a] = Na * user->Function(x,y);
+    F[a] = Na * f;
   }
   return 0;
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal x = p->point[0];
+  PetscReal y = p->point[1];
+  PetscScalar f = user->Function(x,y);
+
+  PetscScalar u;
+  IGAPointFormValue(p,U,&u);
+
+  PetscReal e = PetscAbsScalar(u - f);
+  S[0] = e*e;
+
+  return 0;
+}
+
 #undef __FUNCT__
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
   AppCtx user;
 
-  PetscInt choice=0;
-  const char *choicelist[] = {"paraboloid", "peaks", 0};
-  PetscInt N[2] = {16,16}, nN = 2; 
-  PetscInt p[2] = { 2, 2}, np = 2;
-  PetscInt C[2] = {-1,-1}, nC = 2;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Projection2D Options","IGA");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 = PetscOptionsEList("-function","2D function",__FILE__,choicelist,2,choicelist[choice],&choice,PETSC_NULL);CHKERRQ(ierr);
+  PetscBool flg;
+  PetscBool w[2] = {PETSC_FALSE,PETSC_FALSE}; PetscInt nw = 2;
+  PetscInt  N[2] = {16,16}, nN = 2;
+  PetscInt  p[2] = { 2, 2}, np = 2;
+  PetscInt  C[2] = {-1,-1}, nC = 2;
+  PetscInt  choice=0;
+  const char *choicelist[] = {"paraboloid", "peaks", "hill", "sine", 0};
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","L2 Projection 2D Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBoolArray("-periodic", "periodicity",    __FILE__,w,&nw,&flg);CHKERRQ(ierr);
+  if (flg && nw==0) w[nw++] = PETSC_TRUE;
+  ierr = PetscOptionsIntArray ("-N","number of elements",__FILE__,N,&nN,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray ("-p","polynomial order",  __FILE__,p,&np,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray ("-C","continuity order",  __FILE__,C,&nC,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEList    ("-function","2D function",__FILE__,choicelist,4,choicelist[choice],&choice,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  if (nw == 1) w[1] = w[0];
   if (nN == 1) N[1] = N[0];
   if (np == 1) p[1] = p[0];
   if (nC == 1) C[1] = C[0];
   switch (choice) {
   case 0: user.Function = Paraboloid; break;
   case 1: user.Function = Peaks;      break;
+  case 2: user.Function = Hill;       break;
+  case 3: user.Function = Sine;       break;
   }
 
   IGA iga;
   for (i=0; i<2; i++) {
     IGAAxis axis;
     ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
-  ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis,N[i],-1.0,1.0,C[i]);CHKERRQ(ierr);
+    ierr = IGAAxisSetPeriodic(axis,w[i]);CHKERRQ(ierr);
+    ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],-1.0,1.0,C[i]);CHKERRQ(ierr);
   }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
-  ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
-  
+  PetscScalar error = 0;
+  ierr = IGAFormScalar(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);
+  if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);}
+
+  PetscBool draw = PETSC_FALSE;
+  ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
+  if (draw) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);

demo/L2Projection.c

 #define __FUNCT__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen;
-  IGAPointGetSizes(p,0,&nen,0);
+  PetscInt nen = p->nen;
 
   PetscReal x[3] = {0,0,0};
   IGAPointFormPoint(p,x);
   PetscScalar f = Function(x[0],x[1],x[2]);
 
-  const PetscReal *N;
-  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N);
+  const PetscReal *N = (typeof(N)) p->shape[0];
 
   PetscInt a,b;
   for (a=0; a<nen; a++) {
   PetscInt C[3] = {-1,-1,-1};
   PetscInt n1=3, n2=3, n3=3;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","L2Projection Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,NULL);CHKERRQ(ierr);
   n1 = n2 = n3 = dim;
-  ierr = PetscOptionsIntArray("-N","number of elements",     __FILE__,N,&n1,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray("-p","polynomial order",       __FILE__,p,&n2,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsIntArray("-C","global continuity order",__FILE__,C,&n3,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-N","number of elements",     __FILE__,N,&n1,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-p","polynomial order",       __FILE__,p,&n2,NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-C","global continuity order",__FILE__,C,&n3,NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
   if (n1<3) N[2] = N[0]; if (n1<2) N[1] = N[0];
   if (n2<3) p[2] = p[0]; if (n2<2) p[1] = p[0];
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  ierr = IGASetUserSystem(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,PETSC_NULL);CHKERRQ(ierr);
+  ierr = IGAFormScalar(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);
+/*
+  This code solves the Laplace problem in one of the following ways:
+
+  1) On the parametric unit domain [0,1]^dim (default)
+
+    To solve on the parametric domain, do not specify a geometry
+    file. You may change the discretization by altering the dimension
+    of the space (-iga_dim), the number of uniform elements in each
+    direction (-iga_elements), the polynomial order (-iga_degree),
+    and the continuity (-iga_continuity).
+
+  2) On a geometry
+
+    If a geometry file is specified (-iga_geometry), the
+    discretization will be what is read in from the geometry and is
+    not editable from the command line.
+
+  Note that the boundary conditions for this problem are such that
+  the solution is always u(x)=1 (unit Dirichlet on the left side and
+  free Neumann on the right). The error in the solution may be
+  computed by using the -print_error command.
+
+*/
+
 #include "petiga.h"
 
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
 #undef  __FUNCT__
-#define __FUNCT__ "SystemLaplace"
-PetscErrorCode SystemLaplace(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+#define __FUNCT__ "SystemGalerkin"
+PetscErrorCode SystemGalerkin(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+  const PetscReal (*N1)[dim] = (typeof(N1)) p->shape[1];
 
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
-
-  PetscInt a,b,i;
+  PetscInt a,b;
   for (a=0; a<nen; a++) {
-    for (b=0; b<nen; b++) {
-      PetscScalar Kab = 0.0;
-      for (i=0; i<dim; i++)
-        Kab += N1[a*dim+i]*N1[b*dim+i];
-      K[a*nen+b] = Kab;
-    }
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
+    F[a] = 0.0;
   }
   return 0;
 }
 
-#undef  __FUNCT__
-#define __FUNCT__ "SystemPoisson"
-PetscErrorCode SystemPoisson(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+PETSC_STATIC_INLINE
+PetscReal DEL2(PetscInt dim,const PetscReal a[dim][dim])
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
-
-  const PetscReal *N;
-  IGAPointGetShapeFuns(p,0,&N);
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
-
-  PetscInt a,b,i;
-  for (a=0; a<nen; a++) {
-    for (b=0; b<nen; b++) {
-      PetscScalar Kab = 0.0;
-      for (i=0; i<dim; i++)
-        Kab += N1[a*dim+i]*N1[b*dim+i];
-      K[a*nen+b] = Kab;
-    }
-    F[a] = 1.0*N[a];
-  }
-  return 0;
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i][i];
+  return s;
 }
 
 #undef  __FUNCT__
 #define __FUNCT__ "SystemCollocation"
 PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+  const PetscReal (*N2)[dim][dim] = (typeof(N2)) p->shape[2];
 
-  PetscInt Nb[3] = {0,0,0};
-  Nb[0] = p->parent->parent->axis[0]->nnp;
-  Nb[1] = p->parent->parent->axis[1]->nnp;
-  Nb[2] = p->parent->parent->axis[2]->nnp;
-
-  const PetscReal *N0,(*N1)[dim],(*N2)[dim][dim];
-  IGAPointGetBasisFuns(p,0,(const PetscReal**)&N0);
-  IGAPointGetBasisFuns(p,1,(const PetscReal**)&N1);
-  IGAPointGetBasisFuns(p,2,(const PetscReal**)&N2);
-  
-  PetscInt a,i;
-  PetscReal n[3] = {0,0,0};
-  PetscBool Dirichlet=PETSC_FALSE,Neumann=PETSC_FALSE;
-  for (i=0; i<dim; i++) if (p->parent->ID[i] == 0) Dirichlet=PETSC_TRUE;
-  if (!Dirichlet) { 
-    for (i=0; i<dim; i++) {
-      if (p->parent->ID[i] == Nb[i]-1) {
-	Neumann=PETSC_TRUE;
-	n[i] = 1;
-	i=dim;
-      }
-    }
-  }
-
-  if(Dirichlet){
-    for (a=0; a<nen; a++) K[a] += N0[a];
-    F[0] = 1.0;
-  }else if(Neumann){
-    for (a=0; a<nen; a++) 
-      for (i=0; i<dim; i++) 
-	K[a] += N1[a][i]*n[i];    
-    F[0] = 0.0;
-  }else{
-    for (a=0; a<nen; a++) 
-      for (i=0; i<dim; i++) 
-	K[a] += -N2[a][i][i];
-  }
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] += -DEL2(dim,N2[a]);
+  F[0] = 0.0;
   return 0;
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "ErrorLaplace"
-PetscErrorCode ErrorLaplace(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
 {
   PetscScalar u;
   IGAPointFormValue(p,U,&u);
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
-  /*
-    This code solves the Laplace problem in one of the following ways:
-    
-    1) On the parametric unit domain [0,1]^dim (default)
-    
-      To solve on the parametric domain, do not specify a geometry
-      file. You may change the discretization by altering the
-      dimension of the space (-dim), the number of uniform elements in
-      each direction (-iga_elements), the polynomial order