Commits

Lisandro Dalcin committed 3ba4e3d

More work in collocation and update demos

Comments (0)

Files changed (6)

demo/BoundaryIntegral.c

 /*
   This code solves the Laplace problem where the boundary conditions
-  can be changed from Nuemann to Dirichlet via the commandline. While
+  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.
   return s;
 }
 
-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__ "System"
 PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
   return 0;
 }
 
+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)
   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, PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt ("-side","Neuman BC side",     __FILE__,user.side,&user.side,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-print_error","Prints the L2 error of the solution",__FILE__,print_error,&print_error,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-check_error","Checks the L2 error of the solution",__FILE__,check_error,&check_error,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,PETSC_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,3);CHKERRQ(ierr);}
+  if (iga->dim < 1) {ierr = IGASetDim(iga,2);CHKERRQ(ierr);}
 
   PetscInt dim;
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
-
-  IGABoundary bnd;
   PetscInt d =  !user.side;
   PetscInt n = !!user.side;
   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,PETSC_NULL);CHKERRQ(ierr);
   } else {
+    IGABoundary bnd;
     PetscInt dir,side;
     for (dir=0; dir<dim; dir++) {
       for (side=0; side<2; side++) {
   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);
 
-  PetscBool draw = PETSC_TRUE;
-  if (draw && dim <= 2) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);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);
   return s;
 }
 
-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__ "SystemGalerkin"
 PetscErrorCode SystemGalerkin(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
   return 0;
 }
 
+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)
   for (a=0; a<nen; a++)
     K[a] += -DEL2(dim,N2[a]);
   F[0] = 0.0;
-
   return 0;
 }
 
 
   // Setup options
 
-  PetscInt  i;
   PetscBool collocation = PETSC_FALSE;
   PetscBool print_error = PETSC_FALSE;
+  PetscBool check_error = PETSC_FALSE;
   PetscBool save = PETSC_FALSE;
   PetscBool draw = PETSC_FALSE;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Laplace Options","IGA");CHKERRQ(ierr);
   ierr = PetscOptionsBool("-collocation","Enable to use collocation",__FILE__,collocation,&collocation,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-print_error","Prints the L2 error of the solution",__FILE__,print_error,&print_error,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-check_error","Checks the L2 error of the solution",__FILE__,check_error,&check_error,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-save","Save the solution to file",__FILE__,save,&save,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+  if (collocation) {ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr);}
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
-  if (iga->dim < 1) {ierr = IGASetDim(iga,3);CHKERRQ(ierr);}
+  if (iga->dim < 1) {ierr = IGASetDim(iga,2);CHKERRQ(ierr);}
   ierr = IGASetUp(iga);CHKERRQ(ierr);
-  if (collocation) {ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr);}
 
   // Set boundary conditions
 
-  PetscInt dim;
+  PetscInt dim,dir;
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
-  for (i=0; i<dim; i++) {
+  for (dir=0; dir<dim; dir++) {
     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 = IGABoundarySetLoad(bnd,0,0.0);CHKERRQ(ierr);
+    PetscScalar value = 1.0;
+    PetscScalar load  = 0.0;
+    ierr = IGAGetBoundary(iga,dir,0,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
+    ierr = IGAGetBoundary(iga,dir,1,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetLoad(bnd,0,load);CHKERRQ(ierr);
   }
 
   // Assemble
   if (!iga->collocation) {
     ierr = IGASetUserSystem(iga,SystemGalerkin,PETSC_NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
+    ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
   } else {
     ierr = IGASetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_FALSE);CHKERRQ(ierr);
 
   // Various post-processing options
 
-  if (print_error) {
-    PetscScalar error = 0;
-    ierr = IGAFormScalar(iga,x,1,&error,Error,PETSC_NULL);CHKERRQ(ierr);
-    error = PetscSqrtReal(PetscRealPart(error));
-    ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);
-  }
-  if (save) {
-    ierr = IGAWrite(iga,"LaplaceGeometry.dat");CHKERRQ(ierr);
-    ierr = IGAWriteVec(iga,x,"LaplaceSolution.dat");CHKERRQ(ierr);
-  }
-  if (draw && dim <= 2) {
-    ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
-  }
+  PetscScalar error = 0;
+  ierr = IGAFormScalar(iga,x,1,&error,Error,PETSC_NULL);CHKERRQ(ierr);
+  error = PetscSqrtReal(PetscRealPart(error));
+
+  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);}
+  if (save)        {ierr = IGAWrite(iga,"LaplaceGeometry.dat");CHKERRQ(ierr);}
+  if (save)        {ierr = IGAWriteVec(iga,x,"LaplaceSolution.dat");CHKERRQ(ierr);}
 
   // Cleanup
 
 
 #define pi M_PI
 
-
 PetscReal Exact(PetscReal x[3])
 {
   return sin(pi*x[0]) + sin(pi*x[1]) + sin(pi*x[2]);
   return pi*pi * Exact(x);
 }
 
+PetscReal Flux(PetscInt dir,PetscInt side)
+{
+  return -pi;
+}
+
+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)
+#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;
 
   PetscReal x[3] = {0,0,0};
   IGAPointFormPoint(p,x);
   PetscReal f = Forcing(x);
 
-  const PetscReal *N0,*N1;
-  IGAPointGetShapeFuns(p,0,&N0);
-  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;
+  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] = N0[a]*f;
   }
   return 0;
 }
 
+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;
+
+  PetscReal x[3] = {0,0,0};
+  IGAPointFormPoint(p,x);
+  PetscReal f = Forcing(x);
+
+  const PetscReal (*N2)[dim][dim];
+  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] += -DEL2(dim,N2[a]);
+  F[0] = f;
+  return 0;
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "Error"
 PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
   PetscErrorCode ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
-  PetscInt dim = 3;
+  PetscBool print_error = PETSC_FALSE;
+  PetscBool check_error = PETSC_FALSE;
+  PetscBool draw = PETSC_FALSE;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Neumann Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-dim","dimension",__FILE__,dim,&dim,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-print_error","Prints the L2 error of the solution",__FILE__,print_error,&print_error,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-check_error","Checks the L2 error of the solution",__FILE__,check_error,&check_error,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
-  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-
-  IGABoundary bnd;
-  PetscInt dir,side;
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  if (iga->dim < 1) {ierr = IGASetDim(iga,2);CHKERRQ(ierr);}
+  PetscInt dim,dir,side;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   for (dir=0; dir<dim; dir++) {
     for (side=0; side<2; side++) {
-      PetscScalar load = -pi;
+      IGABoundary bnd;
+      PetscScalar load = Flux(dir,side);
       ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
       ierr = IGABoundarySetLoad(bnd,0,load);CHKERRQ(ierr);
     }
   }
-
-  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   Mat A;
   Vec x,b;
+  KSP ksp;
   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 = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+  if (!iga->collocation) {
+    ierr = IGASetUserSystem(iga,SystemGalerkin,PETSC_NULL);CHKERRQ(ierr);
+    ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
+    ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
+    ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
+  } else {
+    ierr = IGASetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr);
+  }
 
-  KSP ksp;
-  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
-  ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
-  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
-  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   {
     MPI_Comm comm;
     MatNullSpace nsp;
     ierr = KSPSetNullSpace(ksp,nsp);CHKERRQ(ierr);
     ierr = MatNullSpaceDestroy(&nsp);CHKERRQ(ierr);
   }
+
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
   {
-    PetscReal vmin;
+    PetscReal vmin; /* this is a hack */
     ierr = VecMin(x,0,&vmin);CHKERRQ(ierr);
     ierr = VecShift(x,-vmin);CHKERRQ(ierr);
   }
 
-  PetscScalar se = 0;
-  ierr = IGAFormScalar(iga,x,1,&se,Error,PETSC_NULL);CHKERRQ(ierr);
-  PetscReal e = PetscSqrtReal(PetscRealPart(se));
+  PetscScalar error;
+  ierr = IGAFormScalar(iga,x,1,&error,Error,PETSC_NULL);CHKERRQ(ierr);
+  error = PetscSqrtReal(PetscRealPart(error));
 
-  PetscBool error = PETSC_FALSE;
-  ierr = PetscOptionsGetBool(0,"-error",&error,0);CHKERRQ(ierr);
-  if (error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"Error=%G\n",e);CHKERRQ(ierr);}
-  else if (e>1e-4) SETERRQ1(PETSC_COMM_WORLD,1,"Error=%G\n",e);
-
-  if (dim < 3) {
-    PetscBool draw = PETSC_FALSE;
-    ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
-    if (draw) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
-  }
+  if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"Error=%G\n",error);CHKERRQ(ierr);}
+  if (check_error) {if (error>1e-4) SETERRQ1(PETSC_COMM_WORLD,1,"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);
 #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__ "System"
-PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+#define __FUNCT__ "SystemGalerkin"
+PetscErrorCode SystemGalerkin(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt a,b,i,nen=p->nen,dim = p->dim;
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+
   const PetscReal *N0,(*N1)[dim];
   IGAPointGetShapeFuns(p,0,(const PetscReal**)&N0);
   IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+
+  PetscInt a,b;
   for (a=0; a<nen; a++) {
-    for (b=0; b<nen; b++) {
-      for (i=0;i<dim;i++) {
-	K[a*nen+b] += N1[a][i]*N1[b][i];
-      }
-    }
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
     F[a] = N0[a] * 1.0;
   }
   return 0;
 }
 
+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];
+  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] += -DEL2(dim,N2[a]);
+  F[0] = 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);
 
+  PetscBool draw = PETSC_FALSE;
+  PetscBool save = PETSC_FALSE;
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Laplace Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsBool("-save","Save the solution to file",                        __FILE__,save,&save,PETSC_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);}
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
-  IGABoundary bnd;
-  PetscInt    dim,dir,side;
-  PetscScalar value = 1.0;
+  PetscInt dim,dir,side;
   ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
   for (dir=0; dir<dim; dir++) {
     for (side=0; side<2; side++) {
+      IGABoundary bnd;
+      PetscScalar value = 1.0;
       ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
       ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
     }
   ierr = 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,SystemGalerkin,PETSC_NULL);CHKERRQ(ierr);
+    ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
+    ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE);CHKERRQ(ierr);
+  } else {
+    ierr = IGASetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr);
+    ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_FALSE);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);
 
+  if (draw&&dim<3) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+  if (save)        {ierr = IGAWrite(iga,"PoissonGeometry.dat");CHKERRQ(ierr);}
+  if (save)        {ierr = IGAWriteVec(iga,x,"PoissonSolution.dat");CHKERRQ(ierr);}
+
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);
 runex6d_9:
 	-@${MPIEXEC} -n 9 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
 runex7a_1:
-	-@${MPIEXEC} -n 1 ./Neumann ${OPTS} -dim 1
+	-@${MPIEXEC} -n 1 ./Neumann ${OPTS} -iga_dim 1
 runex7a_4:
-	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -dim 1
+	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -iga_dim 1
 runex7b_1:
-	-@${MPIEXEC} -n 1 ./Neumann ${OPTS} -dim 2
+	-@${MPIEXEC} -n 1 ./Neumann ${OPTS} -iga_dim 2
 runex7b_4:
-	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -dim 2
+	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -iga_dim 2
 runex7c_4:
-	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -dim 3
+	-@${MPIEXEC} -n 4 ./Neumann ${OPTS} -iga_dim 3
 runex8_1:
 	-@${MPIEXEC} -n 1 ./ElasticRod ${OPTS} -ts_max_steps 10
 runex8_4:
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidLogicalCollectiveBool(iga,collocation,2);
+  if (iga->setupstage > 0 && iga->collocation != collocation)
+    SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_ARG_WRONGSTATE,
+            "Cannot change collocation after IGASetUp()");
+  iga->collocation = collocation;
   if (collocation && iga->setup) {
     PetscInt i, dim = iga->dim;
     PetscBool periodic = PETSC_FALSE;
     iga->node_iterator->collocation = PETSC_TRUE;
     ierr = IGAElementInit(iga->node_iterator,iga);CHKERRQ(ierr);
   }
-  iga->collocation = collocation;
   PetscFunctionReturn(0);
 }
 
   {
     PetscBool flg;
     PetscInt  i,nw,nl;
+    const char *prefix = 0;
+    PetscBool collocation = iga->collocation;
     IGABasisType btype[3] = {IGA_BASIS_BSPLINE,IGA_BASIS_BSPLINE,IGA_BASIS_BSPLINE};
     PetscBool    wraps[3] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE };
     PetscInt  np,procs[3] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
     PetscInt  dof = (iga->dof > 0) ? iga->dof : 1;
     PetscInt  order = iga->order;
 
+    ierr = IGAGetOptionsPrefix(iga,&prefix);CHKERRQ(ierr);
+
     for (i=0; i<dim; i++) {
       procs[i] = iga->proc_sizes[i];
       wraps[i] = iga->axis[i]->periodic;
     if (flg) {ierr = IGASetDof(iga,dof);CHKERRQ(ierr);}
     dof = (iga->dof > 0) ? iga->dof : 1;
 
+    /* Collocation */
+    ierr = PetscOptionsBool("-iga_collocation","Use collocation","IGASetUseCollocation",collocation,&collocation,&flg);CHKERRQ(ierr);
+    if (flg) {ierr = IGASetUseCollocation(iga,collocation);CHKERRQ(ierr);}
+    if (iga->collocation) {ierr = IGAOptionsReject(prefix,"-iga_quadrature");CHKERRQ(ierr);}
+
     /* Processor grid */
     ierr = PetscOptionsIntArray("-iga_processors","Processor grid","IGASetProcessors",procs,(np=dim,&np),&flg);CHKERRQ(ierr);
     if (flg) for (i=0; i<np; i++) {
     /* Geometry */
     ierr = PetscOptionsString("-iga_geometry","Specify IGA geometry file","IGARead",filename,filename,sizeof(filename),&flg);CHKERRQ(ierr);
     if (flg) { /* load from file */
-      const char *prefix = 0;
-      ierr = IGAGetOptionsPrefix(iga,&prefix);CHKERRQ(ierr);
       ierr = IGAOptionsReject(prefix,"-iga_elements"  );CHKERRQ(ierr);
       ierr = IGAOptionsReject(prefix,"-iga_degree"    );CHKERRQ(ierr);
       ierr = IGAOptionsReject(prefix,"-iga_continuity");CHKERRQ(ierr);
 
   setupcalled:
 
-    /* Collocation */ {
-      PetscBool collocation = iga->collocation;
-      ierr = PetscOptionsBool("-iga_collocation","Use collocation","IGASetCollocation",collocation,&collocation,&flg);CHKERRQ(ierr);
-      if (flg) {ierr = IGASetUseCollocation(iga,collocation);CHKERRQ(ierr);}
-    }
-
     /* Matrix and Vector type */
     if (iga->dof == 1) {ierr = PetscStrcpy(mtype,MATAIJ);CHKERRQ(ierr);}
     if (iga->vectype)  {ierr = PetscStrncpy(vtype,iga->vectype,sizeof(vtype));CHKERRQ(ierr);}
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.