Commits

Nathan Collier committed 86840ab

cleanup pass at demos

  • Participants
  • Parent commits 8af4ba4

Comments (0)

Files changed (13)

File 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;

File 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 {
   IGAPointGetSizes(p,0,&nen,0);
   IGAPointGetDims(p,&dim,0,0);
 
-  const PetscReal *N;
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,0,&N);
-  IGAPointGetShapeFuns(p,1,&N1);
+  const PetscReal *N0,(*N1)[dim];
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
 
   PetscInt a,b,i;
   for (a=0; a<nen; a++) {
       PetscScalar diffusion = 0.0;
       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];
+        diffusion += N1[a][i]*N1[b][i];
+	advection += user->adv[i]*N1[b][i];
       }
-      advection *= N[a];
+      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,PETSC_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;

File demo/BoundaryIntegral.c

+/*
+  This code solves the Laplace problem where the boundary conditions
+  can be changed from Nuemann 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"
 
 #undef  __FUNCT__
   return 0;
 }
 
-
-
 #undef  __FUNCT__
 #define __FUNCT__ "Neumann"
 PetscErrorCode Neumann(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)

File demo/Bratu.c

+/* 
+   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 { 

File demo/ClassicalShell.c

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

File 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];

File 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

File demo/LoggChallenge.c

+/*
+This code was written in response to a challenge posed by Anders Logg
+on Google+. The challenge:
+
+   Solve the partial diferential equation
+
+      -Laplacian(u) = f
+
+   with homogeneous Dirichlet boundary conditions on the unit square for
+
+      f(x,y) = 2 pi^2 sin(pi*x) * sin(pi*y).
+
+   Who can compute a solution with an (L2) error smaller than 10^-6?
+*/
 #include "petiga.h"
 
-/*
-
-Solve the partial diferential equation
-
-   -Laplacian(u) = f
-
-with homogeneous Dirichlet boundary conditions on the unit square for
-
-   f(x,y) = 2 pi^2 sin(pi*x) * sin(pi*y).
-
-+ Who can obtain the smallest error?
-+ Who can compute a solution with an error smaller than 10^-6?
-
-*/
-
 PetscReal Forcing(PetscReal x, PetscReal y)
 {
   PetscReal pi = M_PI;

File demo/PatternFormation.c

+/*
+  This code solves the following system of PDEs
+  
+     u_t = D_1 \nabla^2 u + f(u,v)
+     v_t = D_2 \nabla^2 v + g(u,v)
+
+     f(u,v) = alpha*u*(1-tau1*v**2) + v*(1-tau2*u);
+     g(u,v) = beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);
+
+  and highlights how to use implicit/explicit capabilities of PetIGA.
+*/
 #include "petiga.h"
 
-/*
-
-u_t = D_1 \nabla^2 u + f(u,v)
-v_t = D_2 \nabla^2 v + g(u,v)
-
-f(u,v) = alpha*u*(1-tau1*v**2) + v*(1-tau2*u);
-g(u,v) = beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);
-
-*/
-
 #define EXPLICIT 1
 typedef struct {
   PetscBool IMPLICIT;

File demo/Poisson.c

 #include "petiga.h"
 
-PetscScalar Forcing(const PetscReal x[3],void *ctx)
-{
-  return 1.0;
-}
-
 #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);
-
-  PetscReal x[3] = {0,0,0};
-  IGAPointFormPoint(p,x);
-
-  const PetscReal *N0,*N1;
-  IGAPointGetShapeFuns(p,0,&N0);
-  IGAPointGetShapeFuns(p,1,&N1);
-
-  PetscInt a,b,i;
+  PetscInt a,b,i,nen=p->nen,dim = p->dim;
+  const PetscReal *N0,(*N1)[dim];
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
   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 (i=0;i<dim;i++) {
+	K[a*nen+b] += N1[a][i]*N1[b][i];
+      }
     }
-    F[a] = N0[a]*Forcing(x,ctx);
+    F[a] = N0[a] * 1.0;
   }
   return 0;
 }
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
-  PetscErrorCode ierr;
+  PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
-  PetscInt i,j;
-  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;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson 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 = 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;
-
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
-  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
-  ierr = IGASetDof(iga,dof);CHKERRQ(ierr);
-  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);
-    for (j=0; j<2; j++) {
-      IGABoundary bnd;
-      PetscInt    field = 0;
-      PetscScalar value = 0.0;
-      ierr = IGAGetBoundary(iga,i,j,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,field,value);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  IGABoundary bnd;
+  PetscInt    dim,dir,side;
+  PetscScalar value = 1.0;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  for (dir=0; dir<dim; dir++) {
+    for (side=0; side<2; side++) {
+      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
     }
   }
-  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-  ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   Mat A;
   Vec x,b;
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
   ierr = IGASetUserSystem(iga,System,PETSC_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);
 
-  PetscBool draw = PETSC_FALSE;
-  ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
-  if (draw && dim <= 2) {
-    ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
-  }
-
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);

File demo/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];
+  const PetscReal *N0,(*N1)[1];
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na   = N0[a];
-    PetscReal Na_x = N1[a];
+    PetscReal Na_x = N1[a][0];
     for (b=0; b<nen; b++) {
-      PetscReal Nb_x = N1[b];
+      PetscReal Nb_x = N1[b][0];
       K[a*nen+b] = Na_x * Nb_x;
     }
     F[a] = Na * 1.0;
   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 = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis,N,-1.0,+1.0,C);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 = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
+  IGABoundary bnd;
+  PetscInt dir=0,side;
+  PetscScalar value = 1.0;
+  for (side=0; side<2; side++) {
+    ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);    
+  }
+
   Mat A;
   Vec x,b;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);

File demo/Poisson2D.c

 #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];
+  const PetscReal *N0,(*N1)[2];
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na   = N0[a];
   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
-  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,"","Poisson2D 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 = PetscOptionsEnd();CHKERRQ(ierr);
-  if (nN == 1) N[1] = N[0];
-  if (np == 1) p[1] = p[0];
-  if (nC == 1) C[1] = C[0];
-  if (C[0] == -1) C[0] = p[0]-1;
-  if (C[1] == -1) C[1] = p[1]-1;
-
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,2);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-
-  PetscInt i;
-  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 = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+  
   IGABoundary bnd;
   PetscInt dir,side;
   PetscScalar value = 1.0;
     }
   }
 
-  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-  ierr = IGASetUp(iga);CHKERRQ(ierr);
-
   Mat A;
   Vec x,b;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);

File demo/Poisson3D.c

 #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];
+  const PetscReal *N0,(*N1)[3];
+  IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na   = N0[a];
   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
-  PetscInt N[3] = {16,16,16}, nN = 3; 
-  PetscInt p[3] = { 2, 2, 2}, np = 3;
-  PetscInt C[3] = {-1,-1,-1}, nC = 3;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson2D 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 = PetscOptionsEnd();CHKERRQ(ierr);
-  if (nN == 1) N[2] = N[1] = N[0]; if (nN == 2) N[2] = N[0];
-  if (np == 1) p[2] = p[1] = p[0]; if (np == 2) p[2] = p[0];
-  if (nC == 1) C[2] = C[1] = C[0]; if (nC == 2) C[2] = C[0];
-  if (C[0] == -1) C[0] = p[0]-1;
-  if (C[1] == -1) C[1] = p[1]-1;
-  if (C[2] == -1) C[2] = p[2]-1;
-
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,3);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-  
-  PetscInt i;
-  for (i=0; i<3; 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 = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   IGABoundary bnd;
   PetscInt dir,side;
     }
   }
 
-  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-  ierr = IGASetUp(iga);CHKERRQ(ierr);
-
   Mat A;
   Vec x,b;
   ierr = IGACreateMat(iga,&A);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);