Commits

Lisandro Dalcin committed f3d7eb8

Initial import

  • Participants

Comments (0)

Files changed (48)

+===============
+LICENSE: PetIGA
+===============
+
+:Author:       Lisandro Dalcin, Nathaniel Collier
+:Contact:      dalcinl@gmail.com, nathaniel.collier@gmail.com
+:Organization: CONICET_, KAUST_
+
+.. _CONICET:  http://www.conicet.gov.ar/
+.. _KAUST:    http://www.kaust.edu.sa/
+==============
+README: PetIGA
+==============
+
+:Author:       Lisandro Dalcin, Nathaniel Collier
+:Contact:      dalcinl@gmail.com, nathaniel.collier@gmail.com
+:Organization: CONICET_, KAUST_
+
+.. _CONICET:  http://www.conicet.gov.ar/
+.. _KAUST:    http://www.kaust.edu.sa/
+
+
+- To build and install this package, you must meet the following
+  requirements.
+
+  + PETSc_ 3.2
+
+
+.. _PETSc:  http://www.mcs.anl.gov/petsc/

File conf/petigarules

+# -*- makefile -*-
+include ${PETSC_DIR}/conf/rules
+
+chkpetiga_dir:
+	@mypwd=`pwd`; cd ${PETIGA_DIR} 2>&1 > /dev/null; true_PETIGA_DIR=`pwd`; cd $${mypwd} 2>&1 >/dev/null; \
+        newpwd=`echo $${mypwd} | sed "s+$${true_PETIGA_DIR}+DUMMY+g"`;\
+        haspetiga=`echo $${mypwd} | sed "s+petiga-+DUMMY+g"`;\
+        if [ $${mypwd} = $${newpwd} -a $${haspetiga} != $${mypwd} ]; then \
+          echo "**********************Warning*************************" ; \
+          echo "Your PETIGA_DIR may not match the directory you are in";\
+          echo "PETIGA_DIR " $${true_PETIGA_DIR} "Current directory" $${mypwd};\
+          echo "******************************************************" ; \
+        fi

File conf/petigatest

+# -*- makefile -*-
+include ${PETSC_DIR}/conf/test

File conf/petigavariables

+# -*- makefile -*-
+include ${PETSC_DIR}/conf/variables
+
+PETIGA_INCLUDE = -I${PETIGA_DIR}/${PETSC_ARCH}/include -I${PETIGA_DIR}/include
+PETIGA_LIB_DIR = ${PETIGA_DIR}/${PETSC_ARCH}/lib
+PETIGA_LIB     = ${CC_LINKER_SLFLAG}${PETIGA_LIB_DIR} -L${PETIGA_LIB_DIR} -lpetiga ${PETSC_TS_LIB}
+
+CC_INCLUDES = ${PETSC_INCLUDE} ${PETIGA_INCLUDE}
+FC_INCLUDES = ${PETSC_INCLUDE} ${PETIGA_INCLUDE}
+CCPPFLAGS   = ${PETSC_CCPPFLAGS} ${PETIGA_INCLUDE}
+FCPPFLAGS   = ${PETSC_FCPPFLAGS} ${PETIGA_INCLUDE}
+
+INSTALL_LIB_DIR	= ${PETIGA_LIB_DIR}
+LIBNAME  = ${INSTALL_LIB_DIR}/${LIBBASE}.${AR_LIB_SUFFIX}
+SHLIBS   = libpetiga

File demo/CahnHilliard2D.c

+#include "petiga.h"
+
+typedef struct { 
+  PetscReal theta,cbar,alpha;
+  PetscReal L0,lambda;
+} AppCtx;
+
+#undef  __FUNCT__
+#define __FUNCT__ "Mobility"
+void Mobility(AppCtx *user,PetscReal c,PetscReal *M,PetscReal *dM,PetscReal *d2M)
+{
+  if (M)   *M   = c*(1-c);
+  if (dM)  *dM  = 1-2*c;
+  if (d2M) *d2M = -2;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "ChemicalPotential"
+void ChemicalPotential(AppCtx *user,PetscReal c,PetscReal *mu,PetscReal *dmu,PetscReal *d2mu)
+{
+  if (mu) {
+    (*mu)  = 0.5/user->theta*log(c/(1-c))+1-2*c;
+    (*mu) *= user->L0*user->L0/user->lambda;
+  }
+  if (dmu) {
+    (*dmu)  = 0.5/user->theta*1.0/(c*(1-c)) - 2;
+    (*dmu) *= user->L0*user->L0/user->lambda;
+  }
+  if (d2mu) {
+    (*d2mu)  = -0.5/user->theta*(1-2*c)/(c*c*(1-c)*(1-c));
+    (*d2mu) *= user->L0*user->L0/user->lambda;
+  }
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Residual"
+PetscErrorCode Residual(IGAPoint p,PetscReal dt,PetscReal shift,
+                        PetscReal t,const PetscScalar *U,const PetscScalar *V,
+                        PetscScalar *R,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+
+  PetscScalar c,c_t;
+  IGAPointInterpolate(p,0,U,&c);
+  IGAPointInterpolate(p,0,V,&c_t);
+
+  PetscReal M,dM;
+  Mobility(user,c,&M,&dM,NULL);
+  PetscReal dmu;
+  ChemicalPotential(user,c,NULL,&dmu,NULL);
+
+  PetscScalar c1[2],c2[2][2];
+  IGAPointInterpolate(p,1,U,&c1[0]);
+  IGAPointInterpolate(p,2,U,&c2[0][0]);
+  PetscScalar c_x  = c1[0],    c_y  = c1[1];
+  PetscScalar c_xx = c2[0][0], c_yy = c2[1][1];
+
+  PetscReal *N0 = p->shape[0];
+  PetscReal (*N1)[2] = (PetscReal (*)[2]) p->shape[1];
+  PetscReal (*N2)[2][2] = (PetscReal (*)[2][2]) p->shape[2];
+
+  PetscInt a,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na    = N0[a];
+    PetscReal Na_x  = N1[a][0];
+    PetscReal Na_y  = N1[a][1];
+    PetscReal Na_xx = N2[a][0][0];
+    PetscReal Na_yy = N2[a][1][1];
+    /* ----- */
+    PetscScalar Ra  = 0;
+    // Na * c_t
+    Ra += Na * c_t; 
+    // grad(Na) . ((M*dmu + dM*del2(c))) grad(C)
+    PetscScalar t1 = M*dmu + dM*(c_xx+c_yy);
+    Ra += Na_x * t1 * c_x;
+    Ra += Na_y * t1 * c_y;
+    // del2(Na) * M * del2(c)
+    Ra += (Na_xx+Na_yy) * M * (c_xx+c_yy);
+    /* ----- */
+    R[a] = Ra;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Tangent"
+PetscErrorCode Tangent(IGAPoint p,PetscReal dt,PetscReal shift,
+                       PetscReal t,const PetscScalar *U,const PetscScalar *V,
+                       PetscScalar *K,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+
+  PetscScalar c,c_t;
+  IGAPointInterpolate(p,0,U,&c);
+  IGAPointInterpolate(p,0,V,&c_t);
+
+  PetscReal M,dM,d2M;
+  Mobility(user,c,&M,&dM,&d2M);
+  PetscReal dmu,d2mu;
+  ChemicalPotential(user,c,NULL,&dmu,&d2mu);
+
+  PetscScalar c1[2],c2[2][2];
+  IGAPointInterpolate(p,1,U,&c1[0]);
+  IGAPointInterpolate(p,2,U,&c2[0][0]);
+  PetscScalar c_x  = c1[0],    c_y  = c1[1];
+  PetscScalar c_xx = c2[0][0], c_yy = c2[1][1];
+
+  PetscReal *N0 = p->shape[0];
+  PetscReal (*N1)[2] = (PetscReal (*)[2]) p->shape[1];
+  PetscReal (*N2)[2][2] = (PetscReal (*)[2][2]) p->shape[2];
+
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na    = N0[a];
+    PetscReal Na_x  = N1[a][0];
+    PetscReal Na_y  = N1[a][1];
+    PetscReal Na_xx = N2[a][0][0];
+    PetscReal Na_yy = N2[a][1][1];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb    = N0[b];
+      PetscReal Nb_x  = N1[b][0];
+      PetscReal Nb_y  = N1[b][1];
+      PetscReal Nb_xx = N2[b][0][0];
+      PetscReal Nb_yy = N2[b][1][1];
+      /* ----- */
+      PetscScalar Kab = 0;
+      // shift*Na*Nb
+      Kab += shift*Na*Nb;
+      // grad(Na) . (M*dmu+dM*del2(c)) grad(Nb)
+      PetscScalar t1 = M*dmu + dM*(c_xx+c_yy);
+      Kab += Na_x * t1 * Nb_x;
+      Kab += Na_y * t1 * Nb_y;
+      // grad(Na) . ((dM*dmu+M*d2mu+d2M*del2(c))*Nb + dM*del2(Nb)) grad(C)
+      PetscScalar t2 = (dM*dmu+M*d2mu+d2M*(c_xx+c_yy))*Nb + dM*(Nb_xx+Nb_yy);
+      Kab += Na_x * t2 * c_x;
+      Kab += Na_y * t2 * c_y;
+      // del2(Na) * ((dM*del2(c)*Nb + M*del2(Nb))
+      Kab += (Na_xx+Na_yy) * (dM*(c_xx+c_yy)*Nb + M*(Nb_xx+Nb_yy));
+      /* ----- */
+      K[a*nen+b] = Kab;
+    }
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FormInitialCondition"
+PetscErrorCode FormInitialCondition(AppCtx *user,IGA iga,const char datafile[],Vec C)
+{
+  
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  if (datafile[0] != 0) { /* initial condition from datafile */
+    MPI_Comm comm;
+    PetscViewer viewer;
+    ierr = PetscObjectGetComm((PetscObject)C,&comm);CHKERRQ(ierr);
+    ierr = PetscViewerBinaryOpen(comm,datafile,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
+    ierr = VecLoad(C,viewer);CHKERRQ(ierr);
+    ierr = PetscViewerDestroy(&viewer);
+  } else { /* initial condition is random */
+    PetscRandom rctx;    
+    ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
+    ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
+    ierr = PetscRandomSetInterval(rctx,user->cbar-0.05,user->cbar+0.05);CHKERRQ(ierr); 
+    ierr = PetscRandomSeed(rctx);CHKERRQ(ierr);
+    ierr = VecSetRandom(C,rctx);CHKERRQ(ierr); 
+    ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr); 
+  }
+  PetscFunctionReturn(0); 
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "WriteSolution"
+PetscErrorCode WriteSolution(Vec C, const char pattern[],int number)
+{
+  PetscFunctionBegin;
+  PetscErrorCode  ierr;
+  MPI_Comm        comm;
+  char            filename[256];
+  PetscViewer     viewer;
+
+  PetscFunctionBegin;
+  sprintf(filename,pattern,number);
+  ierr = PetscObjectGetComm((PetscObject)C,&comm);CHKERRQ(ierr);
+  ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
+  ierr = VecView(C,viewer);CHKERRQ(ierr);
+  ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "OutputMonitor"
+PetscErrorCode OutputMonitor(TS ts,PetscInt step,PetscReal t,Vec U,void *mctx)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = WriteSolution(U,"ch%d.dat",step);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  /* Define simulation specific parameters */
+  AppCtx user;
+  user.cbar  = 0.63;   /* average concentration */
+  user.alpha = 3000.0; /* thickess interface parameter */
+  user.theta = 1.5;    /* temperature/critical temperature */
+  user.L0    = 1.0;    /* length scale */
+
+  /* Set discretization options */
+  PetscInt N=64, p=2, C=PETSC_DECIDE;
+  PetscBool output = PETSC_FALSE; 
+  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 = PetscOptionsEnd();CHKERRQ(ierr);
+  if (C == PETSC_DECIDE) C = p-1;
+
+  user.lambda = 1.0/N/N; /* mesh size parameter */
+  
+  if (p < 2 || C < 1) /* Problem requires a p>=2 C1 basis */
+    SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,
+            "Problem requires minimum of p = 2 and C = 1");
+  if (p <= C)         /* Check C < p */
+    SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,
+            "Discretization inconsistent: polynomial order must be greater than degree of continuity");
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,2);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+
+  IGAAxis axis0;
+  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
+  ierr = IGAAxisSetPeriodic(axis0,PETSC_TRUE);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,p,C,N,0.0,1.0);CHKERRQ(ierr);
+  IGAAxis axis1;
+  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
+  ierr = IGAAxisCopy(axis0,axis1);CHKERRQ(ierr);
+
+  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetUserIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
+
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  TS ts;
+  ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
+  ierr = TSSetDuration(ts,10000,1.0);CHKERRQ(ierr);
+  ierr = TSSetTimeStep(ts,1e-10);CHKERRQ(ierr);
+
+  ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
+  ierr = TSAlphaSetRadius(ts,0.5);CHKERRQ(ierr);
+  ierr = TSAlphaSetAdapt(ts,TSAlphaAdaptDefault,PETSC_NULL);CHKERRQ(ierr); 
+
+  if (output) {
+    ierr = TSMonitorSet(ts,OutputMonitor,&user,PETSC_NULL);CHKERRQ(ierr);
+  }
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+
+  PetscReal t; Vec U;
+  ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
+  ierr = FormInitialCondition(&user,iga,initial,U);CHKERRQ(ierr);
+  ierr = TSSolve(ts,U,&t);CHKERRQ(ierr);
+
+  ierr = VecDestroy(&U);CHKERRQ(ierr);
+  ierr = TSDestroy(&ts);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File demo/L2Proj1D.c

+#include "petiga.h"
+
+PetscScalar Line(PetscReal x)
+{
+  return x;
+}
+
+PetscScalar Parabola(PetscReal x)
+{
+  return x*x;
+}
+
+PetscScalar Poly3(PetscReal x)
+{
+  return x*(x-1)*(x+1);
+}
+
+PetscScalar Poly4(PetscReal x)
+{
+  PetscReal a = 1/3.0;
+  return (x-1)*(x-a)*(x+a)*(x+1);
+}
+
+
+typedef struct {
+  PetscScalar (*Function)(PetscReal x);
+} AppCtx;
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal x = p->position[0];
+  PetscReal *N = p->shape[0];
+  
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na = N[a];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb = N[b];
+      K[a*nen+b] = Na * Nb;
+    }
+    F[a] = Na * user->Function(x);
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  AppCtx user;
+
+  PetscInt choice=2;
+  const char *choicelist[] = {"line", "parabola", "poly3", "poly4", 0};
+  PetscInt N=16, p=2, C=PETSC_DECIDE;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Projection1D 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 = PetscOptionsEList("-function","1D function",__FILE__,choicelist,4,choicelist[choice],&choice,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  if (C == PETSC_DECIDE) C = p-1;
+  switch (choice) {
+  case 0: user.Function = Line;     break;
+  case 1: user.Function = Parabola; break;
+  case 2: user.Function = Poly3;    break;
+  case 3: user.Function = Poly4;    break;
+  }
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,1);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  IGAAxis axis0;
+  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,p,C,N,-1.0,1.0);CHKERRQ(ierr);
+
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGAFormSystem(iga,A,b,System,&user);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);
+  
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = MatDestroy(&A);CHKERRQ(ierr);
+  ierr = VecDestroy(&x);CHKERRQ(ierr);
+  ierr = VecDestroy(&b);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File demo/L2Proj2D.c

+#include "petiga.h"
+
+PetscScalar Paraboloid(PetscReal x, PetscReal y)
+{
+  return 1.0 - (x*x + y*y);
+}
+
+PetscScalar Peaks(PetscReal x, PetscReal y)
+{
+  PetscReal X = x*3;
+  PetscReal Y = y*3;
+  return   3 * pow(1-X,2) * exp(-pow(X,2) - pow(Y+1,2))
+    /**/ - 10 * (X/5 - pow(X,3) - pow(Y,5)) * exp(-pow(X,2) - pow(Y,2))
+    /**/ - 1.0/3 * exp(-pow(X+1,2) - pow(Y,2));
+}
+
+typedef struct {
+  PetscScalar (*Function)(PetscReal x, PetscReal y);
+} AppCtx;
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal x = p->position[0];
+  PetscReal y = p->position[1];
+  PetscReal *N = p->shape[0];
+
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na = N[a];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb = N[b];
+      K[a*nen+b] = Na * Nb;
+    }
+    F[a] = Na * user->Function(x,y);
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  AppCtx user;
+
+  PetscInt choice=0;
+  const char *choicelist[] = {"paraboloid", "peaks", 0};
+  PetscInt N=16, p=2, C=PETSC_DECIDE;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Projection2D 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 = PetscOptionsEList("-function","2D function",__FILE__,choicelist,2,choicelist[choice],&choice,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  if (C == PETSC_DECIDE) C = p-1;
+  switch (choice) {
+  case 0: user.Function = Paraboloid; break;
+  case 1: user.Function = Peaks;      break;
+  }
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,2);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  IGAAxis axis0;
+  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,p,C,N,-1.0,1.0);CHKERRQ(ierr);
+  IGAAxis axis1;
+  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
+  ierr = IGAAxisCopy(axis0,axis1);CHKERRQ(ierr);
+
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGAFormSystem(iga,A,b,System,&user);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);
+  
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = MatDestroy(&A);CHKERRQ(ierr);
+  ierr = VecDestroy(&x);CHKERRQ(ierr);
+  ierr = VecDestroy(&b);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File demo/NavierStokesVMS.c

+#include "petiga.h"
+
+typedef struct {
+  PetscReal nu;
+  PetscScalar fx,fy,fz;
+} AppCtx;
+
+#undef __FUNCT__
+#define __FUNCT__ "Tau"
+void Tau(PetscReal GradMap[][3],
+         PetscReal dt,PetscScalar u[],PetscReal nu,
+         PetscScalar *tauM,PetscScalar *tauC)
+{
+  PetscReal C_I = 1.0/12.0;
+
+  PetscReal J[3][3] = {{0}};
+  J[0][0] = 1.0/GradMap[0][0];
+  J[1][1] = 1.0/GradMap[1][1];
+  J[2][2] = 1.0/GradMap[2][2];
+
+  PetscInt i,j,k;
+
+  PetscReal G[3][3] = {{0}};
+  for (i=0;i<3;i++)
+  for (j=0;j<3;j++)
+  for (k=0;k<3;k++)
+    G[i][j] += J[i][k]*J[j][k];
+  PetscReal g[3] = {0};
+  for (i=0;i<3;i++)
+  for (j=0;j<3;j++)
+    g[i] += J[i][j];
+
+  PetscReal G_G = 0;
+  for (i=0;i<3;i++)
+  for (j=0;j<3;j++)
+    G_G += G[i][j]*G[i][j];
+  PetscReal g_g = 0;
+  for (i=0;i<3;i++)
+    g_g += g[i]*g[i];
+
+  PetscScalar u_G_u = 0;
+  for (i=0;i<3;i++)
+  for (j=0;j<3;j++)
+    u_G_u += u[i]*G[i][j]*u[j];
+
+  // Eqn 63
+  *tauM = 4/(dt*dt) + u_G_u + C_I * nu*nu * G_G;
+  *tauM = 1/sqrt(*tauM);
+  // Eqn 64
+  *tauC = (*tauM) * g_g;
+  *tauC = 1/(*tauC);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FineScale"
+void FineScale(const AppCtx *user,PetscScalar tauM,PetscScalar tauC,
+               PetscScalar ux,
+               PetscScalar ux_t,
+               PetscScalar ux_x ,PetscScalar ux_y, PetscScalar ux_z,
+               PetscScalar ux_xx,PetscScalar ux_yy,PetscScalar ux_zz,
+               PetscScalar uy,
+               PetscScalar uy_t,
+               PetscScalar uy_x ,PetscScalar uy_y, PetscScalar uy_z,
+               PetscScalar uy_xx,PetscScalar uy_yy,PetscScalar uy_zz,
+               PetscScalar uz,
+               PetscScalar uz_t,
+               PetscScalar uz_x ,PetscScalar uz_y, PetscScalar uz_z,
+               PetscScalar uz_xx,PetscScalar uz_yy,PetscScalar uz_zz,
+               PetscScalar p,PetscScalar p_x,PetscScalar p_y,PetscScalar p_z,
+               PetscScalar *ux_s,PetscScalar *uy_s,PetscScalar *uz_s,PetscScalar *p_s)
+{
+  // Eqn 61
+  (*ux_s) = ux_t + (ux*ux_x + uy*ux_y + uz*ux_z) + p_x - user->nu*(ux_xx + ux_yy + ux_zz) - user->fx;
+  (*uy_s) = uy_t + (ux*uy_x + uy*uy_y + uz*uy_z) + p_y - user->nu*(uy_xx + uy_yy + uy_zz) - user->fy;
+  (*uz_s) = uz_t + (ux*uz_x + uy*uz_y + uz*uz_z) + p_z - user->nu*(uz_xx + uz_yy + uz_zz) - user->fz;
+  // Eqn 62
+  (*p_s)  = ux_x + uy_y + uz_z;
+  // Eqn 58
+  (*ux_s) *= -tauM;
+  (*uy_s) *= -tauM;
+  (*uz_s) *= -tauM;
+  (*p_s)  *= -tauC;
+}
+
+
+#undef  __FUNCT__
+#define __FUNCT__ "Residual"
+PetscErrorCode Residual(IGAPoint pnt,PetscReal dt,PetscReal shift,
+                        PetscReal t,const PetscScalar *U,const PetscScalar *V,
+                        PetscScalar *Re,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal nu = user->nu;
+
+  PetscScalar u[4],u_t[4];
+  PetscScalar grad_u[4][3];
+  PetscScalar der2_u[4][3][3];
+  IGAPointInterpolate(pnt,0,V,&u_t[0]);
+  IGAPointInterpolate(pnt,0,U,&u[0]);
+  IGAPointInterpolate(pnt,1,U,&grad_u[0][0]);
+  IGAPointInterpolate(pnt,2,U,&der2_u[0][0][0]);
+
+  PetscScalar ux=u[0],ux_t=u_t[0];
+  PetscScalar uy=u[1],uy_t=u_t[1];
+  PetscScalar uz=u[2],uz_t=u_t[2];
+  PetscScalar p =u[3];
+
+  PetscScalar ux_x=grad_u[0][0],ux_y=grad_u[0][1],ux_z=grad_u[0][2];
+  PetscScalar uy_x=grad_u[1][0],uy_y=grad_u[1][1],uy_z=grad_u[1][2];
+  PetscScalar uz_x=grad_u[2][0],uz_y=grad_u[2][1],uz_z=grad_u[2][2];
+  PetscScalar p_x =grad_u[3][0],p_y =grad_u[3][1],p_z =grad_u[3][2];
+
+  PetscScalar ux_xx=der2_u[0][0][0],ux_yy=der2_u[0][1][1],ux_zz=der2_u[0][2][2];
+  PetscScalar uy_xx=der2_u[1][0][0],uy_yy=der2_u[1][1][1],uy_zz=der2_u[1][2][2];
+  PetscScalar uz_xx=der2_u[2][0][0],uz_yy=der2_u[2][1][1],uz_zz=der2_u[2][2][2];
+
+  PetscReal (*GradMap)[3] = (PetscReal(*)[3]) pnt->jacobian;
+  PetscScalar tauM,tauC;
+  Tau(GradMap,dt,u,nu,&tauM,&tauC);
+  PetscScalar ux_s,uy_s,uz_s,p_s;
+  FineScale(user,tauM,tauC,
+            ux,ux_t,ux_x,ux_y,ux_z,ux_xx,ux_yy,ux_zz,
+            uy,uy_t,uy_x,uy_y,uy_z,uy_xx,uy_yy,uy_zz,
+            uz,uz_t,uz_x,uz_y,uz_z,uz_xx,uz_yy,uz_zz,
+            p,p_x,p_y,p_z,
+            &ux_s,&uy_s,&uz_s,&p_s);
+
+  PetscReal  *N0 = pnt->shape[0];
+  PetscReal (*N1)[3] = (PetscReal (*)[3]) pnt->shape[1];
+
+  PetscScalar (*R)[4] = (PetscScalar (*)[4])Re;
+  PetscInt a,nen=pnt->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na    = N0[a];
+    PetscReal Na_x  = N1[a][0];
+    PetscReal Na_y  = N1[a][1];
+    PetscReal Na_z  = N1[a][2];
+    /* ----- */
+    PetscScalar Rux,Ruy,Ruz,Rp;
+    // -L(W)
+    Rux = -Na*user->fx;
+    Ruy = -Na*user->fy;
+    Ruz = -Na*user->fz;
+    Rp  = 0.0;
+    // += B_1(W,U)
+    Rux += Na*ux_t - Na_x*p + nu*( Na_x*( ux_x + ux_x ) + Na_y*( ux_y + uy_x ) + Na_z*( ux_z + uz_x ) );
+    Ruy += Na*uy_t - Na_y*p + nu*( Na_x*( uy_x + ux_y ) + Na_y*( uy_y + uy_y ) + Na_z*( uy_z + uz_y ) );
+    Ruz += Na*uz_t - Na_z*p + nu*( Na_x*( uz_x + ux_z ) + Na_y*( uz_y + uy_z ) + Na_z*( uz_z + uz_z ) );
+    Rp  += Na*( ux_x + uy_y + uz_z );
+    // += Btilde_1(W,U')
+    Rux += - ( Na_x*p_s );
+    Ruy += - ( Na_y*p_s );
+    Ruz += - ( Na_z*p_s );
+    Rp  += - ( Na_x*ux_s + Na_y*uy_s + Na_z*uz_s );
+    // += B_2(W,U,U+U') [after integration by parts, div(u+u')=0]
+    Rux += + Na * ( (ux+ux_s)*ux_x + (uy+uy_s)*ux_y + (uz+uz_s)*ux_z );
+    Ruy += + Na * ( (ux+ux_s)*uy_x + (uy+uy_s)*uy_y + (uz+uz_s)*uy_z );
+    Ruz += + Na * ( (ux+ux_s)*uz_x + (uy+uy_s)*uz_y + (uz+uz_s)*uz_z );
+    // += B_2(W,U',U+U')
+    Rux += - ( Na_x*ux_s*(ux+ux_s) + Na_y*ux_s*(uy+uy_s) + Na_z*ux_s*(uz+uz_s) );
+    Ruy += - ( Na_x*uy_s*(ux+ux_s) + Na_y*uy_s*(uy+uy_s) + Na_z*uy_s*(uz+uz_s) );
+    Ruz += - ( Na_x*uz_s*(ux+ux_s) + Na_y*uz_s*(uy+uy_s) + Na_z*uz_s*(uz+uz_s) );
+    /* ----- */
+    R[a][0] = Rux;
+    R[a][1] = Ruy;
+    R[a][2] = Ruz;
+    R[a][3] = Rp;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "Tangent"
+PetscErrorCode Tangent(IGAPoint pnt,PetscReal dt,PetscReal shift,
+                       PetscReal t,const PetscScalar *U,const PetscScalar *V,
+                       PetscScalar *Ke,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal nu = user->nu;
+
+  PetscScalar u[4];
+  IGAPointInterpolate(pnt,0,U,&u[0]);
+  PetscScalar ux=u[0];
+  PetscScalar uy=u[1];
+  PetscScalar uz=u[2];
+
+  PetscReal (*GradMap)[3] = (PetscReal(*)[3]) pnt->jacobian;
+  PetscScalar tauM,tauC;
+  Tau(GradMap,dt,u,nu,&tauM,&tauC);
+
+  PetscReal *N0 = pnt->shape[0];
+  PetscReal (*N1)[3] = (PetscReal (*)[3]) pnt->shape[1];
+
+  PetscInt a,b,nen=pnt->nen;
+  PetscScalar (*K)[4][nen][4] = (PetscScalar (*)[4][nen][4])Ke;
+  for (a=0; a<nen; a++) {
+    PetscReal Na    = N0[a];
+    PetscReal Na_x  = N1[a][0];
+    PetscReal Na_y  = N1[a][1];
+    PetscReal Na_z  = N1[a][2];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb    = N0[b];
+      PetscReal Nb_x  = N1[b][0];
+      PetscReal Nb_y  = N1[b][1];
+      PetscReal Nb_z  = N1[b][2];
+      /* ----- */
+      PetscInt    i,j;
+      PetscScalar T[4][4];
+      PetscScalar Tii =
+        (+ shift * Na * Nb
+         + Na * (ux * Nb_x + uy * Nb_y + uz * Nb_z)
+         + nu * (Na_x * Nb_x + Na_y * Nb_y + Na_z * Nb_z)
+
+         + tauM * (ux * Na_x + uy * Na_y + uz * Na_z) *
+         /**/     (shift * Nb + (ux * Nb_x + uy * Nb_y + uz * Nb_z))
+         );
+      T[0][0] = /*+ Na * Nb * ux_x*/  +  nu * Na_x * Nb_x  +  tauC * Na_x * Nb_x;
+      T[0][1] = /*+ Na * Nb * ux_y*/  +  nu * Na_y * Nb_x  +  tauC * Na_x * Nb_y;
+      T[0][2] = /*+ Na * Nb * ux_z*/  +  nu * Na_z * Nb_x  +  tauC * Na_x * Nb_z;
+      //
+      T[1][0] = /*+ Na * Nb * uy_x*/  +  nu * Na_x * Nb_y  +  tauC * Na_y * Nb_x;
+      T[1][1] = /*+ Na * Nb * uy_y*/  +  nu * Na_y * Nb_y  +  tauC * Na_y * Nb_y;
+      T[1][2] = /*+ Na * Nb * uy_z*/  +  nu * Na_z * Nb_y  +  tauC * Na_y * Nb_z;
+      //
+      T[2][0] = /*+ Na * Nb * uz_x*/  +  nu * Na_x * Nb_z  +  tauC * Na_z * Nb_x;
+      T[2][1] = /*+ Na * Nb * uz_y*/  +  nu * Na_y * Nb_z  +  tauC * Na_z * Nb_y;
+      T[2][2] = /*+ Na * Nb * uz_z*/  +  nu * Na_z * Nb_z  +  tauC * Na_z * Nb_z;
+      T[0][0] += Tii;
+      T[1][1] += Tii;
+      T[2][2] += Tii;
+      // G as in Eq. (104)
+      T[0][3] = - Na_x * Nb  +  tauM * (ux * Na_x + uy * Na_y + uz * Na_z) * Nb_x;
+      T[1][3] = - Na_y * Nb  +  tauM * (ux * Na_x + uy * Na_y + uz * Na_z) * Nb_y;
+      T[2][3] = - Na_z * Nb  +  tauM * (ux * Na_x + uy * Na_y + uz * Na_z) * Nb_z;
+      // D as in Eq. (106)
+      T[3][0] = + Na * Nb_x  +  tauM * Na_x * (shift * Nb + (ux * Nb_x + uy * Nb_y + uz * Nb_z));
+      T[3][1] = + Na * Nb_y  +  tauM * Na_y * (shift * Nb + (ux * Nb_x + uy * Nb_y + uz * Nb_z));
+      T[3][2] = + Na * Nb_z  +  tauM * Na_z * (shift * Nb + (ux * Nb_x + uy * Nb_y + uz * Nb_z));
+      // L as in Eq. (108)
+      T[3][3] = + tauM * (Na_x * Nb_x + Na_y * Nb_y + Na_z * Nb_z);
+      /* ----- */
+      for (i=0;i<4;i++)
+        for (j=0;j<4;j++)
+          K[a][i][b][j] += T[i][j];
+    }
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FormInitialCondition"
+PetscErrorCode FormInitialCondition(AppCtx *user,IGA iga,const char datafile[],PetscReal t,Vec U)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = VecZeroEntries(U);CHKERRQ(ierr);
+
+  DM da;
+  ierr = IGAGetDofDM(iga,&da);CHKERRQ(ierr);
+  DMDALocalInfo info;
+  ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
+  PetscScalar ****u;
+  ierr = DMDAVecGetArrayDOF(da,U,&u);CHKERRQ(ierr);
+
+  PetscScalar H    = 2;//user->Ly;
+  PetscScalar visc = user->nu;
+  PetscScalar dpdx = user->fx;
+  PetscInt  i,j,k;
+  PetscReal jmax = (info.my-1);
+  for(k=info.zs;k<info.zs+info.zm;k++){
+    for(j=info.ys;j<info.ys+info.ym;j++){
+      for(i=info.xs;i<info.xs+info.xm;i++){
+        PetscReal   y = (j/jmax)*H;
+        PetscScalar ux = 1/(2*visc) * dpdx * y*(H-y);
+        u[k][j][i][0] = ux;
+      }
+    }
+  }
+
+  ierr = DMDAVecRestoreArrayDOF(da,U,&u);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "WriteSolution"
+PetscErrorCode WriteSolution(Vec C, const char pattern[],int number)
+{
+  PetscFunctionBegin;
+  PetscErrorCode  ierr;
+  MPI_Comm        comm;
+  char            filename[256];
+  PetscViewer     viewer;
+
+  PetscFunctionBegin;
+  sprintf(filename,pattern,number);
+  ierr = PetscObjectGetComm((PetscObject)C,&comm);CHKERRQ(ierr);
+  ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
+  ierr = VecView(C,viewer);CHKERRQ(ierr);
+  ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
+  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+#undef __FUNCT__
+#define __FUNCT__ "OutputMonitor"
+PetscErrorCode OutputMonitor(TS ts,PetscInt step,PetscReal t,Vec U,void *mctx)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  ierr = WriteSolution(U,"nsvms%d.dat",step);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  AppCtx user;
+  user.nu = 1.47200e-4;
+  user.fx = 3.37204e-3;
+  user.fy = 0.0;
+  user.fz = 0.0;
+
+  PetscReal PI = 4.0*atan(1.0);
+  PetscReal Lx = 2.0*PI;
+  PetscReal Ly = 2.0;
+  PetscReal Lz = 2.0/3.0*PI;
+
+  PetscInt N=32, p=2, C=PETSC_DECIDE;
+  PetscBool output = PETSC_FALSE;
+  PetscBool monitor = PETSC_FALSE;
+  char initial[PETSC_MAX_PATH_LEN] = {0};
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","NavierStokes 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("-ns_initial","Load initial solution from file",__FILE__,initial,initial,sizeof(initial),PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-ns_output","Enable output files",__FILE__,output,&output,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-ns_monitor","Compute and show statistics of solution",__FILE__,monitor,&monitor,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  if (C == PETSC_DECIDE) C = p-1;
+
+  if (p < 2 || C < 1) /* Problem requires a p>=2 C1 basis */
+    SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,
+            "Problem requires minimum of p = 2 and C = 1");
+  if (p <= C)         /* Check C < p */
+    SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,
+            "Discretization inconsistent: polynomial order must be greater than degree of continuity");
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,3);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,4);CHKERRQ(ierr);
+
+  IGAAxis axis0;
+  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
+  ierr = IGAAxisSetOrder(axis0,p);CHKERRQ(ierr);
+  ierr = IGAAxisSetPeriodic(axis0,PETSC_TRUE);CHKERRQ(ierr);
+  //ierr = IGAAxisInitUniform(axis0,p,C,N,-0.5*Lx,0.5*Lx);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,p,C,N,0,Lx);CHKERRQ(ierr);
+  IGAAxis axis1;
+  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
+  ierr = IGAAxisSetOrder(axis1,p);CHKERRQ(ierr);
+  //ierr = IGAAxisInitUniform(axis1,p,C,N,-0.5*Ly,0.5*Ly);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis1,p,C,N,0,Ly);CHKERRQ(ierr);
+  IGAAxis axis2;
+  ierr = IGAGetAxis(iga,2,&axis2);CHKERRQ(ierr);
+  ierr = IGAAxisSetOrder(axis2,p);CHKERRQ(ierr);
+  ierr = IGAAxisSetPeriodic(axis2,PETSC_TRUE);CHKERRQ(ierr);
+  //ierr = IGAAxisInitUniform(axis2,p,C,N,-0.5*Lz,0.5*Lz);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis2,p,C,N,0,Lz);CHKERRQ(ierr);
+
+  IGABoundary bnd;
+  PetscInt dir=1,side,field;
+  for (side=0;side<2;side++) {
+    ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+    for (field=0;field<3;field++) {
+      ierr = IGABoundarySetValue(bnd,field,0.0);CHKERRQ(ierr);
+    }
+  }
+
+  ierr = IGASetUserIFunction(iga,Residual,&user);CHKERRQ(ierr);
+  ierr = IGASetUserIJacobian(iga,Tangent,&user);CHKERRQ(ierr);
+
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  TS ts;
+  ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
+  ierr = TSSetDuration(ts,1000000,1000.0);CHKERRQ(ierr);
+  ierr = TSSetTimeStep(ts,1.0e-2);CHKERRQ(ierr);
+
+  ierr = TSSetType(ts,TSALPHA);CHKERRQ(ierr);
+  ierr = TSAlphaSetRadius(ts,0.5);CHKERRQ(ierr);
+
+  if (output) {
+    ierr = TSMonitorSet(ts,OutputMonitor,&user,PETSC_NULL);CHKERRQ(ierr);
+  }
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+
+  PetscReal t=0; Vec U;
+  ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
+  ierr = FormInitialCondition(&user,iga,initial,t,U);CHKERRQ(ierr);
+  ierr = TSSolve(ts,U,&t);CHKERRQ(ierr);
+
+  ierr = VecDestroy(&U);CHKERRQ(ierr);
+  ierr = TSDestroy(&ts);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File demo/Poisson1D.c

+#include "petiga.h"
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscReal *N0 = p->shape[0];
+  PetscReal *N1 = p->shape[1];
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na   = N0[a];
+    PetscReal Na_x = N1[a];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb_x = N1[b];
+      K[a*nen+b] = Na_x * Nb_x;
+    }
+    F[a] = Na * 1.0;
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  PetscInt N=16, p=2, C=PETSC_DECIDE;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson1D 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 = PetscOptionsEnd();CHKERRQ(ierr);
+  if (C == PETSC_DECIDE) C = p-1;
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,1);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+
+  IGAAxis axis;
+  ierr = IGAGetAxis(iga,0,&axis);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis,p,C,N,-1.0,+1.0);CHKERRQ(ierr);
+
+  IGABoundary left;
+  IGABoundary right;
+  PetscScalar value = 1.0;
+  ierr = IGAGetBoundary(iga,0,0,&left);CHKERRQ(ierr);
+  ierr = IGABoundarySetValue(left,0,value);CHKERRQ(ierr);
+  ierr = IGAGetBoundary(iga,0,1,&right);CHKERRQ(ierr);
+  ierr = IGABoundarySetValue(right,0,value);CHKERRQ(ierr);
+
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGAFormSystem(iga,A,b,System,PETSC_NULL);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);
+  
+  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
+  ierr = MatDestroy(&A);CHKERRQ(ierr);
+  ierr = VecDestroy(&x);CHKERRQ(ierr);
+  ierr = VecDestroy(&b);CHKERRQ(ierr);
+  ierr = IGADestroy(&iga);CHKERRQ(ierr);
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File demo/Poisson2D.c

+#include "petiga.h"
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscReal *N0 = p->shape[0];
+  PetscReal (*N1)[2] = (PetscReal (*)[2]) p->shape[1];
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na   = N0[a];
+    PetscReal Na_x = N1[a][0];
+    PetscReal Na_y = N1[a][1];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb_x = N1[b][0];
+      PetscReal Nb_y = N1[b][1];
+      K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y;
+    }
+    F[a] = Na * 1.0;
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  PetscInt N=16, p=2, C=PETSC_DECIDE;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson2D 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 = PetscOptionsEnd();CHKERRQ(ierr);
+  if (C == PETSC_DECIDE) C = p-1;
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,2);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+
+  IGAAxis axis0;
+  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,p,C,N,-1.0,1.0);CHKERRQ(ierr);
+  IGAAxis axis1;
+  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
+  ierr = IGAAxisCopy(axis0,axis1);CHKERRQ(ierr);
+
+  IGABoundary bnd;
+  PetscInt dir,side;
+  PetscScalar value = 1.0;
+  for (dir=0; dir<2; dir++) {
+    for (side=0; side<2; side++) {
+      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
+    }
+  }
+
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGAFormSystem(iga,A,b,System,PETSC_NULL);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);
+  
+  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;
+}

File demo/Poisson3D.c

+#include "petiga.h"
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscReal *N0 = p->shape[0];
+  PetscReal (*N1)[3] = (PetscReal (*)[3]) p->shape[1];
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na   = N0[a];
+    PetscReal Na_x = N1[a][0];
+    PetscReal Na_y = N1[a][1];
+    PetscReal Na_z = N1[a][2];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb_x = N1[b][0];
+      PetscReal Nb_y = N1[b][1];
+      PetscReal Nb_z = N1[b][2];
+      K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y + Na_z*Nb_z;
+    }
+    F[a] = Na * 1.0;
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  PetscInt N=16, p=2, C=PETSC_DECIDE;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson3D 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 = PetscOptionsEnd();CHKERRQ(ierr);
+  if (C == PETSC_DECIDE) C = p-1;
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,3);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  
+  IGAAxis axis0;
+  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,p,C,N,-1.0,1.0);CHKERRQ(ierr);
+  IGAAxis axis1;
+  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
+  ierr = IGAAxisCopy(axis0,axis1);CHKERRQ(ierr);
+  IGAAxis axis2;
+  ierr = IGAGetAxis(iga,2,&axis2);CHKERRQ(ierr);
+  ierr = IGAAxisCopy(axis1,axis2);CHKERRQ(ierr);
+
+  IGABoundary bnd;
+  PetscInt dir,side;
+  PetscScalar value = 1.0;
+  for (dir=0; dir<3; dir++) {
+    for (side=0; side<2; side++) {
+      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
+    }
+  }
+
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGAFormSystem(iga,A,b,System,PETSC_NULL);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);
+
+  /*
+  PetscViewer viewer = PETSC_VIEWER_BINARY_WORLD;
+  ierr = VecView(x,viewer);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);
+
+  ierr = PetscFinalize();CHKERRQ(ierr);
+  return 0;
+}

File demo/makefile

+TARGETS = L2Proj1D L2Proj2D \
+          Poisson1D Poisson2D Poisson3D \
+          CahnHilliard2D \
+          NavierStokesVMS
+ALL: ${TARGETS}
+clean::
+	-@${RM} ${TARGETS}
+
+CFLAGS    = #-g3 -Wall -Wextra -Wno-unused-parameter #-Wconversion
+FFLAGS    = #-g3 -Wall -Wextra -fcheck=all
+CPPFLAGS  =
+FPPFLAGS  =
+LOCDIR    = demo
+EXAMPLESC = CahnHilliard2D.c
+EXAMPLESF =
+MANSEC    = IGA
+
+topdir := $(shell cd .. && pwd)
+PETIGA_DIR ?= $(topdir)
+include ${PETIGA_DIR}/conf/petigavariables
+include ${PETIGA_DIR}/conf/petigarules
+
+L2Proj1D: L2Proj1D.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
+L2Proj2D: L2Proj2D.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
+Poisson1D: Poisson1D.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
+Poisson2D: Poisson2D.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
+Poisson3D: Poisson3D.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
+CahnHilliard2D: CahnHilliard2D.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
+NavierStokesVMS: NavierStokesVMS.o chkopts
+	${CLINKER} -o $@ $< ${PETIGA_LIB}
+	${RM} -f $<
+
+runex1a_1:
+	-@${MPIEXEC} -n 1 ./L2Proj1D -nox -malloc_debug -malloc_dump
+runex1a_4:
+	-@${MPIEXEC} -n 4 ./L2Proj1D -nox -malloc_debug -malloc_dump
+runex1b_1:
+	-@${MPIEXEC} -n 1 ./L2Proj2D -nox -malloc_debug -malloc_dump
+runex1b_4:
+	-@${MPIEXEC} -n 4 ./L2Proj2D -nox -malloc_debug -malloc_dump
+
+runex2a_1:
+	-@${MPIEXEC} -n 1 ./Poisson1D -nox -malloc_debug -malloc_dump
+runex2a_4:
+	-@${MPIEXEC} -n 4 ./Poisson1D -nox -malloc_debug -malloc_dump
+runex2b_1:
+	-@${MPIEXEC} -n 1 ./Poisson2D -nox -malloc_debug -malloc_dump
+runex2b_4:
+	-@${MPIEXEC} -n 4 ./Poisson2D -nox -malloc_debug -malloc_dump
+runex2c_1:
+	-@${MPIEXEC} -n 1 ./Poisson3D -nox -malloc_debug -malloc_dump
+runex2c_4:
+	-@${MPIEXEC} -n 4 ./Poisson3D -nox -malloc_debug -malloc_dump
+runex4_1:
+	-@${MPIEXEC} -n 1 ./CahnHilliard2D -ts_max_steps 2 -malloc_debug -malloc_dump
+runex4_4:
+	-@${MPIEXEC} -n 2 ./CahnHilliard2D -ts_max_steps 2 -malloc_debug -malloc_dump
+
+L2Proj1D := L2Proj1D.PETSc runex1a_1 runex1a_4 L2Proj1D.rm
+L2Proj2D := L2Proj2D.PETSc runex1b_1 runex1b_4 L2Proj2D.rm
+Poisson1D := Poisson1D.PETSc runex2a_1 runex2a_4 Poisson1D.rm
+Poisson2D := Poisson2D.PETSc runex2b_1 runex2b_4 Poisson2D.rm
+Poisson3D := Poisson3D.PETSc runex2c_1 runex2c_4 Poisson3D.rm
+CahnHilliard2D := CahnHilliard2D.PETSc runex4_1 runex4_4 CahnHilliard2D.rm
+
+L2Proj       := $(L2Proj1D) $(L2Proj2D) $(L2Proj3D)
+Poisson      := $(Poisson1D) $(Poisson2D) $(Poisson3D)
+CahnHilliard := $(CahnHilliard2D) $(CahnHilliard3D)
+
+TESTEXAMPLES_C := $(L2Proj) $(Poisson) $(CahnHilliard)
+TESTEXAMPLES_F :=
+TESTEXAMPLES_FORTRAN:=$(TESTEXAMPLES_F)
+testexamples:
+	-@${OMAKE} tree ACTION=testexamples_C PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR}
+testfortran:
+	-@if [ "${FC}" != "" ]; then \
+            ${OMAKE} tree ACTION=testexamples_Fortran PETSC_ARCH=${PETSC_ARCH} PETSC_DIR=${PETSC_DIR} PETIGA_DIR=${PETIGA_DIR}; \
+          fi
+
+include ${PETIGA_DIR}/conf/petigatest

File docs/makefile

+DIRS     = 
+LOCDIR   = docs/
+
+include ${PETIGA_DIR}/conf/petigavariables
+include ${PETIGA_DIR}/conf/petigarules
+include ${PETIGA_DIR}/conf/petigatest

File docs/manual/INSTALL.txt

+Install
+=======
+
+XXX To be written ...
+
+.. Local Variables:
+.. mode: rst
+.. End:

File docs/manual/Makefile

+# Makefile for Sphinx documentation
+#
+
+# You can set these variables from the command line.
+SPHINXOPTS    =
+SPHINXBUILD   = sphinx-build
+PAPER         =
+BUILDDIR      = _build
+
+# Internal variables.
+PAPEROPT_a4     = -D latex_paper_size=a4
+PAPEROPT_letter = -D latex_paper_size=letter
+ALLSPHINXOPTS   = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) .
+
+.PHONY: help clean html dirhtml singlehtml pickle json htmlhelp qthelp devhelp epub latex latexpdf text man changes linkcheck doctest
+
+help:
+	@echo "Please use \`make <target>' where <target> is one of"
+	@echo "  html       to make standalone HTML files"
+	@echo "  dirhtml    to make HTML files named index.html in directories"
+	@echo "  singlehtml to make a single large HTML file"
+	@echo "  pickle     to make pickle files"
+	@echo "  json       to make JSON files"
+	@echo "  htmlhelp   to make HTML files and a HTML help project"
+	@echo "  qthelp     to make HTML files and a qthelp project"
+	@echo "  devhelp    to make HTML files and a Devhelp project"
+	@echo "  epub       to make an epub"
+	@echo "  latex      to make LaTeX files, you can set PAPER=a4 or PAPER=letter"
+	@echo "  latexpdf   to make LaTeX files and run them through pdflatex"
+	@echo "  text       to make text files"
+	@echo "  man        to make manual pages"
+	@echo "  changes    to make an overview of all changed/added/deprecated items"
+	@echo "  linkcheck  to check all external links for integrity"
+	@echo "  doctest    to run all doctests embedded in the documentation (if enabled)"
+
+clean:
+	-rm -rf $(BUILDDIR)
+
+html:
+	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
+	@echo
+	@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
+
+dirhtml:
+	$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
+	@echo
+	@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
+
+singlehtml:
+	$(SPHINXBUILD) -b singlehtml $(ALLSPHINXOPTS) $(BUILDDIR)/singlehtml
+	@echo
+	@echo "Build finished. The HTML page is in $(BUILDDIR)/singlehtml."
+
+pickle:
+	$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
+	@echo
+	@echo "Build finished; now you can process the pickle files."
+
+json:
+	$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
+	@echo
+	@echo "Build finished; now you can process the JSON files."
+
+htmlhelp:
+	$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
+	@echo
+	@echo "Build finished; now you can run HTML Help Workshop with the" \
+	      ".hhp project file in $(BUILDDIR)/htmlhelp."
+
+qthelp:
+	$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
+	@echo
+	@echo "Build finished; now you can run "qcollectiongenerator" with the" \
+	      ".qhcp project file in $(BUILDDIR)/qthelp, like this:"
+	@echo "# qcollectiongenerator $(BUILDDIR)/qthelp/PetIGA.qhcp"
+	@echo "To view the help file:"
+	@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/PetIGA.qhc"
+
+devhelp:
+	$(SPHINXBUILD) -b devhelp $(ALLSPHINXOPTS) $(BUILDDIR)/devhelp
+	@echo
+	@echo "Build finished."
+	@echo "To view the help file:"
+	@echo "# mkdir -p $$HOME/.local/share/devhelp/PetIGA"
+	@echo "# ln -s $(BUILDDIR)/devhelp $$HOME/.local/share/devhelp/PetIGA"
+	@echo "# devhelp"
+
+epub:
+	$(SPHINXBUILD) -b epub $(ALLSPHINXOPTS) $(BUILDDIR)/epub
+	@echo
+	@echo "Build finished. The epub file is in $(BUILDDIR)/epub."
+
+latex:
+	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+	@echo
+	@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
+	@echo "Run \`make' in that directory to run these through (pdf)latex" \
+	      "(use \`make latexpdf' here to do that automatically)."
+
+latexpdf:
+	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
+	@echo "Running LaTeX files through pdflatex..."
+	make -C $(BUILDDIR)/latex all-pdf
+	@echo "pdflatex finished; the PDF files are in $(BUILDDIR)/latex."
+
+text:
+	$(SPHINXBUILD) -b text $(ALLSPHINXOPTS) $(BUILDDIR)/text
+	@echo
+	@echo "Build finished. The text files are in $(BUILDDIR)/text."
+
+man:
+	$(SPHINXBUILD) -b man $(ALLSPHINXOPTS) $(BUILDDIR)/man
+	@echo
+	@echo "Build finished. The manual pages are in $(BUILDDIR)/man."
+
+changes:
+	$(SPHINXBUILD) -b changes $(ALLSPHINXOPTS) $(BUILDDIR)/changes
+	@echo
+	@echo "The overview file is in $(BUILDDIR)/changes."
+
+linkcheck:
+	$(SPHINXBUILD) -b linkcheck $(ALLSPHINXOPTS) $(BUILDDIR)/linkcheck
+	@echo
+	@echo "Link check complete; look for any errors in the above output " \
+	      "or in $(BUILDDIR)/linkcheck/output.txt."
+
+doctest:
+	$(SPHINXBUILD) -b doctest $(ALLSPHINXOPTS) $(BUILDDIR)/doctest
+	@echo "Testing of doctests in the sources finished, look at the " \
+	      "results in $(BUILDDIR)/doctest/output.txt."

File docs/manual/OVERVIEW.txt

+Overview
+========
+
+XXX To be written ...
+
+.. Local Variables:
+.. mode: rst
+.. End:

File docs/manual/TUTORIAL.txt

+Tutorial
+========
+
+XXX To be written ...
+
+.. Local Variables:
+.. mode: rst
+.. End:

File docs/manual/abstract.txt

+.. topic:: Abstract
+
+   This document describes PetIGA ...
+
+.. Local Variables:
+.. mode: rst
+.. End:

File docs/manual/conf.py

+# -*- coding: utf-8 -*-
+#
+# PetIGA documentation build configuration file, created by
+# sphinx-quickstart on Fri Sep 23 11:42:32 2011.
+#
+# This file is execfile()d with the current directory set to its containing dir.
+#
+# Note that not all possible configuration values are present in this
+# autogenerated file.
+#
+# All configuration values have a default; values that are commented out
+# serve to show the default.
+
+import sys, os
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#sys.path.insert(0, os.path.abspath('.'))
+
+# -- General configuration -----------------------------------------------------
+
+# If your documentation needs a minimal Sphinx version, state it here.
+#needs_sphinx = '1.0'
+
+# Add any Sphinx extension module names here, as strings. They can be extensions
+# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
+#extensions = ['sphinx.ext.pngmath', 'sphinx.ext.jsmath']
+extensions = ['sphinx.ext.pngmath']
+
+# Add any paths that contain templates here, relative to this directory.
+#templates_path = ['_templates']
+
+# The suffix of source filenames.
+source_suffix = '.txt'
+
+# The encoding of source files.
+#source_encoding = 'utf-8-sig'
+
+# The master toctree document.
+master_doc = 'index'
+
+# General information about the project.
+project = u'PetIGA'
+copyright = u'2011, L. Dalcin, N. Collier'
+
+# The version info for the project you're documenting, acts as replacement for
+# |version| and |release|, also used in various other places throughout the
+# built documents.
+#
+# The short X.Y version.
+version = '0.1'
+# The full version, including alpha/beta/rc tags.
+release = '0.1'
+
+# The language for content autogenerated by Sphinx. Refer to documentation
+# for a list of supported languages.
+#language = None
+
+# There are two options for replacing |today|: either, you set today to some
+# non-false value, then it is used:
+#today = ''
+# Else, today_fmt is used as the format for a strftime call.
+#today_fmt = '%B %d, %Y'
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+exclude_patterns = ['_build']
+
+# The reST default role (used for this markup: `text`) to use for all documents.
+#default_role = None
+
+# If true, '()' will be appended to :func: etc. cross-reference text.
+#add_function_parentheses = True
+
+# If true, the current module name will be prepended to all description
+# unit titles (such as .. function::).
+#add_module_names = True
+
+# If true, sectionauthor and moduleauthor directives will be shown in the
+# output. They are ignored by default.
+#show_authors = False
+
+# The name of the Pygments (syntax highlighting) style to use.
+pygments_style = 'sphinx'
+
+# A list of ignored prefixes for module index sorting.
+#modindex_common_prefix = []
+
+
+# -- Options for HTML output ---------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+html_theme = 'default'
+
+# Theme options are theme-specific and customize the look and feel of a theme
+# further.  For a list of options available for each theme, see the
+# documentation.
+#html_theme_options = {}
+
+# Add any paths that contain custom themes here, relative to this directory.
+#html_theme_path = []
+
+# The name for this set of Sphinx documents.  If None, it defaults to
+# "<project> v<release> documentation".
+#html_title = None
+
+# A shorter title for the navigation bar.  Default is the same as html_title.
+#html_short_title = None
+
+# The name of an image file (relative to this directory) to place at the top
+# of the sidebar.
+#html_logo = None
+
+# The name of an image file (within the static path) to use as favicon of the
+# docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
+# pixels large.
+#html_favicon = None
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+#html_static_path = ['_static']
+
+# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
+# using the given strftime format.
+#html_last_updated_fmt = '%b %d, %Y'
+
+# If true, SmartyPants will be used to convert quotes and dashes to
+# typographically correct entities.
+#html_use_smartypants = True
+
+# Custom sidebar templates, maps document names to template names.
+#html_sidebars = {}
+
+# Additional templates that should be rendered to pages, maps page names to
+# template names.
+#html_additional_pages = {}
+
+# If false, no module index is generated.
+#html_domain_indices = True
+
+# If false, no index is generated.
+#html_use_index = True
+
+# If true, the index is split into individual pages for each letter.
+#html_split_index = False
+
+# If true, links to the reST sources are added to the pages.
+#html_show_sourcelink = True
+
+# If true, "Created using Sphinx" is shown in the HTML footer. Default is True.
+#html_show_sphinx = True
+
+# If true, "(C) Copyright ..." is shown in the HTML footer. Default is True.
+#html_show_copyright = True
+
+# If true, an OpenSearch description file will be output, and all pages will
+# contain a <link> tag referring to it.  The value of this option must be the
+# base URL from which the finished HTML is served.
+#html_use_opensearch = ''
+
+# This is the file name suffix for HTML files (e.g. ".xhtml").
+#html_file_suffix = None
+
+# Output file base name for HTML help builder.
+htmlhelp_basename = 'PetIGAdoc'
+
+
+# -- Options for LaTeX output --------------------------------------------------
+
+# The paper size ('letter' or 'a4').
+#latex_paper_size = 'letter'
+latex_paper_size = 'a4'
+
+# The font size ('10pt', '11pt' or '12pt').
+#latex_font_size = '10pt'
+
+# Grouping the document tree into LaTeX files. List of tuples
+# (source start file, target name, title, author, documentclass [howto/manual]).
+latex_documents = [
+  ('manual', 'PetIGA.tex', u'PetIGA Documentation',
+   u'Lisandro Dalcin, Nathaniel Collier', 'howto'),
+]
+
+# The name of an image file (relative to this directory) to place at the top of
+# the title page.
+#latex_logo = None
+
+# For "manual" documents, if this is true, then toplevel headings are parts,
+# not chapters.
+#latex_use_parts = False
+
+# If true, show page references after internal links.
+#latex_show_pagerefs = False
+
+# If true, show URL addresses after external links.
+#latex_show_urls = False
+
+# Additional stuff for the LaTeX preamble.
+#latex_preamble = ''
+
+# Documents to append as an appendix to all manuals.
+#latex_appendices = []
+
+# If false, no module index is generated.
+#latex_domain_indices = True
+
+latex_additional_files = ['sphinxfix.sty']
+latex_elements = {
+    'printmodindex': '',
+    'printindex': '',
+    'preamble' : r'\usepackage{sphinxfix}',
+    }
+
+# -- Options for manual page output --------------------------------------------
+
+# One entry per manual page. List of tuples
+# (source start file, name, description, authors, manual section).
+man_pages = [
+    ('index', 'petiga', u'PetIGA Documentation',
+     [u'Lisandro Dalcin, Nathaniel Collier'], 1)
+]

File docs/manual/index.txt

+======
+PetIGA
+======
+
+:Authors:      Nathaniel Collier, Lisandro Dalcin
+:Contact:       nathaniel.collier@gmail.com, dalcinl@gmail.com
+:Organization: KAUST_, CONICET_
+:Date:         |today|
+
+.. include:: abstract.txt
+
+Contents
+========
+
+.. include:: toctree.txt
+
+.. include:: links.txt
+
+Index
+=====
+
+* :ref:`genindex`
+* :ref:`search`
+
+.. Local Variables:
+.. mode: rst
+.. End:

File docs/manual/links.txt

+.. _CONICET:  http://www.conicet.gov.ar/
+.. _CIMEC:    http://www.cimec.org.ar/
+.. _KAUST:    http://www.kaust.edu.sa/
+
+.. _MPI:      http://www.mpi-forum.org/
+.. _PETSc:    http://www.mcs.anl.gov/petsc/
+
+.. Local Variables:
+.. mode: rst
+.. End:

File docs/manual/make.bat

+@ECHO OFF
+
+REM Command file for Sphinx documentation
+
+if "%SPHINXBUILD%" == "" (
+	set SPHINXBUILD=sphinx-build
+)
+set BUILDDIR=_build
+set ALLSPHINXOPTS=-d %BUILDDIR%/doctrees %SPHINXOPTS% .
+if NOT "%PAPER%" == "" (
+	set ALLSPHINXOPTS=-D latex_paper_size=%PAPER% %ALLSPHINXOPTS%
+)
+
+if "%1" == "" goto help
+
+if "%1" == "help" (
+	:help
+	echo.Please use `make ^<target^>` where ^<target^> is one of
+	echo.  html       to make standalone HTML files
+	echo.  dirhtml    to make HTML files named index.html in directories
+	echo.  singlehtml to make a single large HTML file
+	echo.  pickle     to make pickle files
+	echo.  json       to make JSON files
+	echo.  htmlhelp   to make HTML files and a HTML help project
+	echo.  qthelp     to make HTML files and a qthelp project
+	echo.  devhelp    to make HTML files and a Devhelp project
+	echo.  epub       to make an epub
+	echo.  latex      to make LaTeX files, you can set PAPER=a4 or PAPER=letter
+	echo.  text       to make text files
+	echo.  man        to make manual pages
+	echo.  changes    to make an overview over all changed/added/deprecated items
+	echo.  linkcheck  to check all external links for integrity
+	echo.  doctest    to run all doctests embedded in the documentation if enabled
+	goto end
+)
+
+if "%1" == "clean" (
+	for /d %%i in (%BUILDDIR%\*) do rmdir /q /s %%i
+	del /q /s %BUILDDIR%\*
+	goto end
+)
+
+if "%1" == "html" (
+	%SPHINXBUILD% -b html %ALLSPHINXOPTS% %BUILDDIR%/html
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished. The HTML pages are in %BUILDDIR%/html.
+	goto end
+)
+
+if "%1" == "dirhtml" (
+	%SPHINXBUILD% -b dirhtml %ALLSPHINXOPTS% %BUILDDIR%/dirhtml
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished. The HTML pages are in %BUILDDIR%/dirhtml.
+	goto end
+)
+
+if "%1" == "singlehtml" (
+	%SPHINXBUILD% -b singlehtml %ALLSPHINXOPTS% %BUILDDIR%/singlehtml
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished. The HTML pages are in %BUILDDIR%/singlehtml.
+	goto end
+)
+
+if "%1" == "pickle" (
+	%SPHINXBUILD% -b pickle %ALLSPHINXOPTS% %BUILDDIR%/pickle
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished; now you can process the pickle files.
+	goto end
+)
+
+if "%1" == "json" (
+	%SPHINXBUILD% -b json %ALLSPHINXOPTS% %BUILDDIR%/json
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished; now you can process the JSON files.
+	goto end
+)
+
+if "%1" == "htmlhelp" (
+	%SPHINXBUILD% -b htmlhelp %ALLSPHINXOPTS% %BUILDDIR%/htmlhelp
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished; now you can run HTML Help Workshop with the ^
+.hhp project file in %BUILDDIR%/htmlhelp.
+	goto end
+)
+
+if "%1" == "qthelp" (
+	%SPHINXBUILD% -b qthelp %ALLSPHINXOPTS% %BUILDDIR%/qthelp
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished; now you can run "qcollectiongenerator" with the ^
+.qhcp project file in %BUILDDIR%/qthelp, like this:
+	echo.^> qcollectiongenerator %BUILDDIR%\qthelp\PetIGA.qhcp
+	echo.To view the help file:
+	echo.^> assistant -collectionFile %BUILDDIR%\qthelp\PetIGA.ghc
+	goto end
+)
+
+if "%1" == "devhelp" (
+	%SPHINXBUILD% -b devhelp %ALLSPHINXOPTS% %BUILDDIR%/devhelp
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished.
+	goto end
+)
+
+if "%1" == "epub" (
+	%SPHINXBUILD% -b epub %ALLSPHINXOPTS% %BUILDDIR%/epub
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished. The epub file is in %BUILDDIR%/epub.
+	goto end
+)
+
+if "%1" == "latex" (
+	%SPHINXBUILD% -b latex %ALLSPHINXOPTS% %BUILDDIR%/latex
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished; the LaTeX files are in %BUILDDIR%/latex.
+	goto end
+)
+
+if "%1" == "text" (
+	%SPHINXBUILD% -b text %ALLSPHINXOPTS% %BUILDDIR%/text
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished. The text files are in %BUILDDIR%/text.
+	goto end
+)
+
+if "%1" == "man" (
+	%SPHINXBUILD% -b man %ALLSPHINXOPTS% %BUILDDIR%/man
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Build finished. The manual pages are in %BUILDDIR%/man.
+	goto end
+)
+
+if "%1" == "changes" (
+	%SPHINXBUILD% -b changes %ALLSPHINXOPTS% %BUILDDIR%/changes
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.The overview file is in %BUILDDIR%/changes.
+	goto end
+)
+
+if "%1" == "linkcheck" (
+	%SPHINXBUILD% -b linkcheck %ALLSPHINXOPTS% %BUILDDIR%/linkcheck
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Link check complete; look for any errors in the above output ^
+or in %BUILDDIR%/linkcheck/output.txt.
+	goto end
+)
+
+if "%1" == "doctest" (
+	%SPHINXBUILD% -b doctest %ALLSPHINXOPTS% %BUILDDIR%/doctest
+	if errorlevel 1 exit /b 1
+	echo.
+	echo.Testing of doctests in the sources finished, look at the ^
+results in %BUILDDIR%/doctest/output.txt.
+	goto end
+)
+
+:end

File docs/manual/manual.txt

+======
+PetIGA
+======
+
+.. include:: abstract.txt
+
+.. include:: toctree.txt
+
+.. include:: links.txt
+
+.. Local Variables:
+.. mode: rst
+.. End:

File docs/manual/sphinxfix.sty

+\setcounter{tocdepth}{2}
+
+\pagenumbering{arabic}
+
+\makeatletter
+\renewcommand{\theindex}{
+  \cleardoublepage
+  \phantomsection
+  \py@OldTheindex
+  \addcontentsline{toc}{section}{\indexname}
+}
+\makeatother
+
+\makeatletter
+\renewcommand{\thebibliography}[1]{
+  \cleardoublepage
+  \phantomsection
+  \py@OldThebibliography{1}
+  \addcontentsline{toc}{section}{\bibname}
+}
+\makeatother
+
+\makeatletter
+\renewcommand{\tableofcontents}{
+  \begingroup
+    \parskip = 0mm
+    \py@OldTableofcontents
+  \endgroup
+  \vfill
+  \rule{\textwidth}{1pt}
+  \newpage
+}
+\makeatother
+

File docs/manual/toctree.txt

+.. toctree::
+   :maxdepth: 2
+
+   OVERVIEW
+   INSTALL
+   TUTORIAL
+
+
+.. Local Variables:
+.. mode: rst
+.. End:

File include/makefile

+CFLAGS   =
+FFLAGS   =
+SOURCEH  = petiga.h
+SOURCEC  =
+SOURCEF  =
+OBJSC    =
+OBJSF    =
+LIBBASE  = libpetiga
+DIRS     = #finclude private
+LOCDIR   = include/
+MANSEC   = IGA
+
+include ${PETIGA_DIR}/conf/petigavariables
+include ${PETIGA_DIR}/conf/petigarules
+include ${PETIGA_DIR}/conf/petigatest

File include/petiga.h

+#if !defined(PETIGA_H)
+#define PETIGA_H
+
+/*
+#include "petscconf.h"
+#undef  PETSC_STATIC_INLINE
+#define PETSC_STATIC_INLINE static __inline
+*/
+
+#include <petsc.h>
+#include <petscdmda.h>
+PETSC_EXTERN_CXX_BEGIN
+
+#ifndef PetscMalloc1
+#define PetscMalloc1(m1,t1,r1) (PetscMalloc((m1)*sizeof(t1),(r1)))
+#endif
+ 
+/* ---------------------------------------------------------------- */
+
+typedef struct _n_IGAAxis     *IGAAxis;
+typedef struct _n_IGARule     *IGARule;
+typedef struct _n_IGABasis    *IGABasis;
+typedef struct _n_IGABoundary *IGABoundary;
+
+typedef struct _n_IGAElement  *IGAElement;
+typedef struct _n_IGAPoint    *IGAPoint;
+
+/* ---------------------------------------------------------------- */