1. Lisandro Dalcin
  2. PetIGA

Commits

Nathan Collier  committed 12bff28 Merge

branch merge

  • Participants
  • Parent commits 86840ab, 43b015a
  • Branches default

Comments (0)

Files changed (23)

File demo/BoundaryIntegral.c

View file
  */
 #include "petiga.h"
 
+typedef struct {
+  PetscInt dir;
+  PetscInt side;
+} AppCtx;
+
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
+{
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
+
+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)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
-
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
-
-  PetscInt a,b,i;
+  PetscInt dim = p->dim;
+  PetscInt nen = p->nen;
+  const PetscReal (*N1)[dim]; //*((PetscReal**)&N1) = p->shape[1];
+  IGAPointGetShapeFuns(p,1,(const PetscReal **)&N1);
+  PetscInt a,b;
   for (a=0; a<nen; a++) {
-    for (b=0; b<nen; b++) {
-      PetscScalar Kab = 0.0;
-      for (i=0; i<dim; i++)
-        Kab += N1[a*dim+i]*N1[b*dim+i];
-      K[a*nen+b] = Kab;
-    }
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
+    F[a] = 0.0;
   }
   return 0;
 }
 #define __FUNCT__ "Neumann"
 PetscErrorCode Neumann(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscReal *N0 = p->shape[0];
-  PetscInt a,nen=p->nen;
-  for (a=0; a<nen; a++) {
-    PetscReal Na   = N0[a];
-    F[a] = Na * 1.0;
-  }
+  PetscInt nen = p->nen;
+  const PetscReal *N0 = p->shape[0];
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    F[a] = N0[a] * 1.0;
   return 0;
 }
 
-typedef struct { 
-  PetscInt dir;
-  PetscInt side;
-} AppCtx;
+#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]; //*((PetscReal**)&N2) = p->shape[2];
+  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
+
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = -DEL2(dim,N2[a]);
+  F[0] = 0.0;
+
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "DirichletCollocation"
+PetscErrorCode DirichletCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen = p->nen;
+  const PetscReal *N0 = p->shape[0];
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = N0[a];
+  F[0] = 1.0;
+  return 0;
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "NeumannCollocation"
+PetscErrorCode NeumannCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
+  const PetscReal (*N1)[dim];
+  const PetscReal *normal = p->normal;
+  IGAPointGetShapeFuns(p,1,(const PetscReal**)&N1);
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] = DOT(dim,N1[a],normal);
+  F[0] = 1.0;
+  return 0;
+}
+PetscErrorCode Neumann0Collocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{(void)NeumannCollocation(p,K,F,ctx); F[0] = 0.0; return 0;}
 
 #undef  __FUNCT__
 #define __FUNCT__ "Error"
   IGAPointFormValue(p,U,&u);
   PetscReal x;
   if (user->side == 0)
-    x = 1 - p->point[user->dir];
+    x = 1 - p->point[user->dir] + 1;
   else
-    x = p->point[user->dir];
+    x = p->point[user->dir] + 1;
   PetscReal e = u - x;
   S[0] = e*e;
   return 0;
   user.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); 
+  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);
   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);}
+
+  PetscInt dim;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
 
   IGABoundary bnd;
-  PetscInt d = !user.side; 
+  PetscInt d =  !user.side;
   PetscInt n = !!user.side;
-  ierr = IGAGetBoundary(iga,user.dir,d,&bnd);CHKERRQ(ierr);
-  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-  ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr); 
-  ierr = IGABoundarySetUserSystem(bnd,Neumann,PETSC_NULL);CHKERRQ(ierr);
-
-  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  if (!iga->collocation) {
+    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 {
+    PetscInt dir,side;
+    for (dir=0; dir<dim; dir++) {
+      for (side=0; side<2; side++) {
+        ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+        ierr = IGABoundarySetUserSystem(bnd,Neumann0Collocation,PETSC_NULL);CHKERRQ(ierr);
+      }
+    }
+    ierr = IGAGetBoundary(iga,user.dir,d,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetUserSystem(bnd,DirichletCollocation,PETSC_NULL);CHKERRQ(ierr);
+    ierr = IGAGetBoundary(iga,user.dir,n,&bnd);CHKERRQ(ierr);
+    ierr = IGABoundarySetUserSystem(bnd,NeumannCollocation,PETSC_NULL);CHKERRQ(ierr);
+  }
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   Mat A;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  if (!iga->collocation) {
+    ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  } else {
+    ierr = IGASetUserSystem(iga,SystemCollocation,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);
 
-  PetscInt dim;
-  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
-  if (dim <= 2) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
-
   PetscScalar error = 0;
   ierr = IGAFormScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
   error = PetscSqrtReal(PetscRealPart(error));
   ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);
-  
+
+  PetscBool draw = PETSC_TRUE;
+  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);
   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/Bratu.c

View file
   if (steady) {
     SNES snes;
     ierr = IGACreateSNES(iga,&snes);CHKERRQ(ierr);
+    ierr = SNESSetTolerances(snes,PETSC_DEFAULT,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
     ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
     ierr = SNESSolve(snes,0,x);CHKERRQ(ierr);
     ierr = SNESDestroy(&snes);CHKERRQ(ierr);
   } else {
     TS ts;
+    SNES snes;
     ierr = IGACreateTS(iga,&ts);CHKERRQ(ierr);
     ierr = TSSetType(ts,TSTHETA);CHKERRQ(ierr);
     ierr = TSSetDuration(ts,10000,0.1);CHKERRQ(ierr);
     ierr = TSSetTimeStep(ts,0.01);CHKERRQ(ierr);
+    ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
+    ierr = SNESSetTolerances(snes,PETSC_DEFAULT,1e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
     ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
 #if PETSC_VERSION_LE(3,3,0)
     ierr = TSSolve(ts,x,PETSC_NULL);CHKERRQ(ierr);

File demo/Laplace.c

View file
+/*
+  This code solves the Laplace problem in one of the following ways:
+
+  1) On the parametric unit domain [0,1]^dim (default)
+
+    To solve on the parametric domain, do not specify a geometry
+    file. You may change the discretization by altering the dimension
+    of the space (-iga_dim), the number of uniform elements in each
+    direction (-iga_elements), the polynomial order (-iga_degree),
+    and the continuity (-iga_continuity).
+
+  2) On a geometry
+
+    If a geometry file is specified (-iga_geometry), the
+    discretization will be what is read in from the geometry and is
+    not editable from the command line.
+
+  Note that the boundary conditions for this problem are such that
+  the solution is always u(x)=1 (unit Dirichlet on the left side and
+  free Neumann on the right). The error in the solution may be
+  computed by using the -print_error command.
+
+*/
+
 #include "petiga.h"
 
-#undef  __FUNCT__
-#define __FUNCT__ "SystemLaplace"
-PetscErrorCode SystemLaplace(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+PETSC_STATIC_INLINE
+PetscReal DOT(PetscInt dim,const PetscReal a[],const PetscReal b[])
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
+  PetscInt i; PetscReal s = 0.0;
+  for (i=0; i<dim; i++) s += a[i]*b[i];
+  return s;
+}
 
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
-
-  PetscInt a,b,i;
-  for (a=0; a<nen; a++) {
-    for (b=0; b<nen; b++) {
-      PetscScalar Kab = 0.0;
-      for (i=0; i<dim; i++)
-        Kab += N1[a*dim+i]*N1[b*dim+i];
-      K[a*nen+b] = Kab;
-    }
-  }
-  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__ "SystemPoisson"
-PetscErrorCode SystemPoisson(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+#define __FUNCT__ "SystemGalerkin"
+PetscErrorCode SystemGalerkin(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
 
-  const PetscReal *N;
-  IGAPointGetShapeFuns(p,0,&N);
-  const PetscReal *N1;
-  IGAPointGetShapeFuns(p,1,&N1);
+  const PetscReal (*N1)[dim];
+  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;
-    }
-    F[a] = 1.0*N[a];
+    for (b=0; b<nen; b++)
+      K[a*nen+b] = DOT(dim,N1[a],N1[b]);
+    F[a] = 0.0;
   }
   return 0;
 }
 #define __FUNCT__ "SystemCollocation"
 PetscErrorCode SystemCollocation(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
 {
-  PetscInt nen,dim;
-  IGAPointGetSizes(p,0,&nen,0);
-  IGAPointGetDims(p,&dim,0,0);
+  PetscInt nen = p->nen;
+  PetscInt dim = p->dim;
 
-  PetscInt Nb[3] = {0,0,0};
-  Nb[0] = p->parent->parent->axis[0]->nnp;
-  Nb[1] = p->parent->parent->axis[1]->nnp;
-  Nb[2] = p->parent->parent->axis[2]->nnp;
+  const PetscReal (*N2)[dim][dim];
+  IGAPointGetShapeFuns(p,2,(const PetscReal**)&N2);
 
-  const PetscReal *N0,(*N1)[dim],(*N2)[dim][dim];
-  IGAPointGetBasisFuns(p,0,(const PetscReal**)&N0);
-  IGAPointGetBasisFuns(p,1,(const PetscReal**)&N1);
-  IGAPointGetBasisFuns(p,2,(const PetscReal**)&N2);
-  
-  PetscInt a,i;
-  PetscReal n[3] = {0,0,0};
-  PetscBool Dirichlet=PETSC_FALSE,Neumann=PETSC_FALSE;
-  for (i=0; i<dim; i++) if (p->parent->ID[i] == 0) Dirichlet=PETSC_TRUE;
-  if (!Dirichlet) { 
-    for (i=0; i<dim; i++) {
-      if (p->parent->ID[i] == Nb[i]-1) {
-	Neumann=PETSC_TRUE;
-	n[i] = 1;
-	i=dim;
-      }
-    }
-  }
+  PetscInt a;
+  for (a=0; a<nen; a++)
+    K[a] += -DEL2(dim,N2[a]);
+  F[0] = 0.0;
 
-  if(Dirichlet){
-    for (a=0; a<nen; a++) K[a] += N0[a];
-    F[0] = 1.0;
-  }else if(Neumann){
-    for (a=0; a<nen; a++) 
-      for (i=0; i<dim; i++) 
-	K[a] += N1[a][i]*n[i];    
-    F[0] = 0.0;
-  }else{
-    for (a=0; a<nen; a++) 
-      for (i=0; i<dim; i++) 
-	K[a] += -N2[a][i][i];
-  }
   return 0;
 }
 
 #undef  __FUNCT__
-#define __FUNCT__ "ErrorLaplace"
-PetscErrorCode ErrorLaplace(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
 {
   PetscScalar u;
   IGAPointFormValue(p,U,&u);
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
-  /*
-    This code solves the Laplace problem in one of the following ways:
-    
-    1) On the parametric unit domain [0,1]^dim (default)
-    
-      To solve on the parametric domain, do not specify a geometry
-      file. You may change the discretization by altering the
-      dimension of the space (-dim), the number of uniform elements in
-      each direction (-iga_elements), the polynomial order
-      (-iga_degree), and the continuity (-iga_continuity). Note that
-      the boundary conditions for this problem are such that the
-      solution is always u(x)=1 (unit Dirichlet on the left side and
-      free Neumann on the right). The error in the solution may be
-      computed by using the -print_error command.
-
-    2) On a geometry
-
-      If a geometry file is specified (-iga_geometry), then the code
-      will solve the Poisson problem on this geometry. The forcing is
-      set to 1 and we have 0 Dirichlet conditions everywhere. We use
-      this mode to test geometries as the solution should display the
-      same symmetry as the geometry. The discretization will be what
-      is read in from the geometry and is not editable from the
-      commandline.
-
-   */
-
   PetscErrorCode ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
   // Setup options
 
-  PetscInt  i; 
-  PetscInt  dim = 3; 
-  PetscBool print_error = PETSC_FALSE; 
-  PetscBool draw = PETSC_FALSE; 
-  PetscBool Collocation = PETSC_FALSE;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Laplace 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("-draw","If dim <= 2, then draw the solution to the screen",__FILE__,draw,&draw,PETSC_NULL);CHKERRQ(ierr); 
-  ierr = PetscOptionsBool("-collocation","Enable to use collocation",__FILE__,Collocation,&Collocation,PETSC_NULL);CHKERRQ(ierr); 
-  ierr = PetscOptionsEnd();CHKERRQ(ierr); 
+  PetscInt  i;
+  PetscBool collocation = PETSC_FALSE;
+  PetscBool print_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("-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);
 
   // Initialize the discretization
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-  ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  if (iga->dim < 1) {ierr = IGASetDim(iga,3);CHKERRQ(ierr);}
   ierr = IGASetUp(iga);CHKERRQ(ierr);
+  if (collocation) {ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr);}
 
   // Set boundary conditions
 
-  if (iga->geometry) {
-    for(i=0; i<dim; i++) {
-      IGABoundary bnd;
-      ierr = IGAGetBoundary(iga,i,0,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-      ierr = IGAGetBoundary(iga,i,1,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-    }
-  }else{
-    for (i=0; i<dim; i++) {
-      IGABoundary bnd;
-      ierr = IGAGetBoundary(iga,i,0,&bnd);CHKERRQ(ierr);
-      ierr = IGABoundarySetValue(bnd,0,1.0);CHKERRQ(ierr);
-    }
+  PetscInt dim;
+  ierr = IGAGetDim(iga,&dim);CHKERRQ(ierr);
+  for (i=0; i<dim; i++) {
+    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);
   }
-  
+
   // Assemble
 
   Mat A;
-  Vec x,b;
+  Vec b,x;
   ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
   ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
-  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
-  if (Collocation){
-    ierr = IGASetUseCollocation(iga,PETSC_TRUE);CHKERRQ(ierr); 
-    ierr = IGASetUserSystem(iga,SystemCollocation,PETSC_NULL);CHKERRQ(ierr); 
-    ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
-  }else{
-    if (iga->geometry){ 
-      ierr = IGASetUserSystem(iga,SystemPoisson,PETSC_NULL);CHKERRQ(ierr); 
-    }else{
-      ierr = IGASetUserSystem(iga,SystemLaplace,PETSC_NULL);CHKERRQ(ierr); 
-    }
-    ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  if (!iga->collocation) {
+    ierr = IGASetUserSystem(iga,SystemGalerkin,PETSC_NULL);CHKERRQ(ierr);
     ierr = MatSetOption(A,MAT_SYMMETRIC,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);
 
   // Solve
 
 
   // Various post-processing options
 
-  if (iga->geometry) {
-    MPI_Comm        comm;
-    PetscViewer     viewer;
-    ierr = PetscObjectGetComm((PetscObject)x,&comm);CHKERRQ(ierr);
-    ierr = PetscViewerBinaryOpen(comm,"solution.dat",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
-    ierr = VecView(x,viewer);CHKERRQ(ierr);
-    ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
-    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
-  }
-
-  if (print_error && !iga->geometry) {
+  if (print_error) {
     PetscScalar error = 0;
-    ierr = IGAFormScalar(iga,x,1,&error,ErrorLaplace,PETSC_NULL);CHKERRQ(ierr);
+    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 (draw && dim <= 2) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);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);
+  }
 
   // Cleanup
 

File demo/NavierStokesKorteweg2D.c

View file
   PetscReal t=0; Vec U;
   ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
   ierr = FormInitialCondition(iga,t,U,&user);CHKERRQ(ierr);
+#if PETSC_VERSION_LE(3,3,0)
   ierr = TSSolve(ts,U,PETSC_NULL);CHKERRQ(ierr);
+#else
+  ierr = TSSolve(ts,U);CHKERRQ(ierr);
+#endif
 
   ierr = VecDestroy(&U);CHKERRQ(ierr);
   ierr = TSDestroy(&ts);CHKERRQ(ierr);

File demo/Richards.c

View file
   ierr = TSAlphaSetRadius(ts,0.5);CHKERRQ(ierr);
   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
 
-  PetscReal t; 
   Vec       U;
   ierr = IGACreateVec(iga,&U);CHKERRQ(ierr);
   ierr = FormInitialCondition(iga,U,&user);CHKERRQ(ierr);
-  ierr = TSSolve(ts,U,&t);CHKERRQ(ierr);
+#if PETSC_VERSION_LE(3,3,0)
+  ierr = TSSolve(ts,U,PETSC_NULL);CHKERRQ(ierr);
+#else
+  ierr = TSSolve(ts,U);CHKERRQ(ierr);
+#endif
 
   ierr = VecDestroy(&U);CHKERRQ(ierr);
   ierr = TSDestroy(&ts);CHKERRQ(ierr);

File demo/makefile

View file
 	  Elasticity3D \
 	  HyperElasticity \
 	  PhaseFieldCrystal2D \
-          BoundaryIntegral \
-          Richards \
-          TwoPhaseTwoComponent \
+      BoundaryIntegral \
+      Richards \
+      TwoPhaseTwoComponent \
 	  ShallowWater \
-          ClassicalShell \
+      ClassicalShell \
 	  LoggChallenge \
-          ElasticRod
+      ElasticRod
 
 ALL: ${TARGETS}
 clean::
 runex6c_4:
 	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady false -ts_max_steps 2
 runex6d_1:
-	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2 -steady -iga_collocation
+	-@${MPIEXEC} -n 1 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
+runex6d_2:
+	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
+runex6d_4:
+	-@${MPIEXEC} -n 4 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
+runex6d_8:
+	-@${MPIEXEC} -n 8 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
+runex6d_9:
+	-@${MPIEXEC} -n 9 ./Bratu ${OPTS} -iga_dim 2 -steady true -iga_collocation
 runex7a_1:
 	-@${MPIEXEC} -n 1 ./Neumann ${OPTS} -dim 1
 runex7a_4:
 Poisson2D := Poisson2D.PETSc runex2b_1 runex2b_4 Poisson2D.rm
 Poisson3D := Poisson3D.PETSc runex2c_1 runex2c_4 Poisson3D.rm
 Neumann := Neumann.PETSc runex7a_1 runex7a_4 runex7b_1 runex7b_4 runex7c_4 Neumann.rm
-Bratu := Bratu.PETSc runex6a_1 runex6a_2 runex6b_1 runex6b_4 runex6c_1 runex6c_4 runex6d_1 Bratu.rm
+Bratu   := \
+Bratu.PETSc \
+runex6a_1 runex6a_2 runex6b_1 runex6b_4 runex6c_1 runex6c_4 \
+runex6d_1 runex6d_2 runex6d_4 runex6d_8 runex6d_9 \
+Bratu.rm
 CahnHilliard2D := CahnHilliard2D.PETSc runex4_1 runex4_4 CahnHilliard2D.rm
 PatternFormation := PatternFormation.PETSc runex5a_1 runex5a_4 runex5b_1 runex5b_4 PatternFormation.rm
 ElasticRod := ElasticRod.PETSc runex8_1 runex8_4 ElasticRod.rm

File include/petiga.h

View file
 PETSC_EXTERN PetscErrorCode IGAAxisGetPeriodic(IGAAxis axis,PetscBool *periodic);
 PETSC_EXTERN PetscErrorCode IGAAxisSetDegree(IGAAxis axis,PetscInt p);
 PETSC_EXTERN PetscErrorCode IGAAxisGetDegree(IGAAxis axis,PetscInt *p);
-PETSC_EXTERN PetscErrorCode IGAAxisSetKnots(IGAAxis axis,PetscInt m,PetscReal U[]);
+PETSC_EXTERN PetscErrorCode IGAAxisSetKnots(IGAAxis axis,PetscInt m,const PetscReal U[]);
 PETSC_EXTERN PetscErrorCode IGAAxisGetKnots(IGAAxis axis,PetscInt *m,PetscReal *U[]);
 PETSC_EXTERN PetscErrorCode IGAAxisGetLimits(IGAAxis axis,PetscReal *Ui,PetscReal *Uf);
 PETSC_EXTERN PetscErrorCode IGAAxisGetSizes(IGAAxis axis,PetscInt *nel,PetscInt *nnp);
+PETSC_EXTERN PetscErrorCode IGAAxisGetSpans(IGAAxis axis,PetscInt *nel,PetscInt *spans[]);
+PETSC_EXTERN PetscErrorCode IGAAxisInit(IGAAxis axis,PetscInt p,PetscInt m,const PetscReal U[]);
 PETSC_EXTERN PetscErrorCode IGAAxisInitBreaks(IGAAxis axis,PetscInt nu,const PetscReal u[],PetscInt C);
 PETSC_EXTERN PetscErrorCode IGAAxisInitUniform(IGAAxis axis,PetscInt N,PetscReal Ui,PetscReal Uf,PetscInt C);
 PETSC_EXTERN PetscErrorCode IGAAxisSetUp(IGAAxis axis);
 PETSC_EXTERN PetscErrorCode IGAPrependOptionsPrefix(IGA iga,const char prefix[]);
 PETSC_EXTERN PetscErrorCode IGAAppendOptionsPrefix(IGA iga,const char prefix[]);
 PETSC_EXTERN PetscErrorCode IGASetFromOptions(IGA iga);
+PETSC_EXTERN PetscErrorCode IGAOptionsAlias(const char name[],const char defval[],const char alias[]);
 
 PETSC_EXTERN PetscErrorCode IGALoad(IGA iga,PetscViewer viewer);
 PETSC_EXTERN PetscErrorCode IGASave(IGA iga,PetscViewer viewer);
   PetscInt  width[3];
   PetscInt  sizes[3];
   PetscInt  ID[3];
+  IGABasis  BD[3];
   PetscBool atboundary;
-  PetscInt  atboundary_id;
+  PetscInt  boundary_id;
   /**/
   PetscInt count;
   PetscInt index;
   PetscInt nsd;
   PetscInt npd;
 
-  IGABasis *BD;
-
   PetscInt    *mapping;   /*[nen]      */
   PetscBool   geometry;
   PetscReal   *geometryX; /*[nen][nsd] */

File src/makefile

View file
 
 #FFLAGS_STD_F03 = -std=f2003
 #FFLAGS_STD_F08 = -std=f2008
+#FFLAGS_STD_F08 = -std=f2008ts
 #FFLAGS_STD = ${FFLAGS_STD_F03}
-#FFLAGS = ${FFLAGS_STD} -pedantic -Wall -Wextra -fcheck=all
+#FFLAGS = ${FFLAGS_STD} -pedantic -Wall -Wextra -Wimplicit-interface -g3 -fcheck=all -fbacktrace
 
 SOURCEH  = ../include/petiga.h petigabl.h petigagrid.h petigapart.h
 SOURCEC  = petiga.c petigareg.c petigaaxis.c petigarule.c petigabasis.c petigabound.c petigaelem.c petigapoint.c petigavec.c petigamat.c petigascl.c petigapcb.c petigapce.c petigapc.c petigaksp.c petigasnes.c petigats.c petigats2.c petigaio.c petigagrid.c petigapart.c snesfdcolor.c tsalpha2.c

File src/petiga.c

View file
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   PetscValidLogicalCollectiveBool(iga,collocation,2);
-  if (collocation && !iga->collocation) {
-    PetscMPIInt size = 1;
-    PetscInt i, dim = (iga->dim > 0) ? iga->dim : 3;
+  if (collocation && iga->setup) {
+    PetscInt i, dim = iga->dim;
     PetscBool periodic = PETSC_FALSE;
-    ierr = MPI_Comm_size(((PetscObject)iga)->comm,&size);CHKERRQ(ierr);
     for (i=0; i<dim; i++) if(iga->axis[i]->periodic) periodic = PETSC_TRUE;
-    if (size > 1) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_SUP,
-                          "Collocation not supported in parallel");
     if (periodic) SETERRQ(((PetscObject)iga)->comm,PETSC_ERR_SUP,
                           "Collocation not supported with periodicity");
   }
-  if (collocation && iga->setup) { /* collocation */
+  if (collocation && iga->setup) {
     PetscInt i;
     for (i=0; i<3; i++) {
       ierr = IGABasisDestroy(&iga->node_basis[i]);CHKERRQ(ierr);
   PetscFunctionReturn(0);
 }
 
+PETSC_STATIC_INLINE
+#undef  __FUNCT__
+#define __FUNCT__ "IGAOptionsReject"
+PetscErrorCode IGAOptionsReject(const char prefix[],const char name[])
+{
+  PetscErrorCode ierr;
+  PetscBool      flag = PETSC_FALSE;
+  PetscFunctionBegin;
+  ierr = PetscOptionsHasName(prefix,name,&flag);CHKERRQ(ierr);
+  if (flag) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Disabled option: %s",name);
+  PetscFunctionReturn(0);
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGASetFromOptions"
 /*@
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   {
     PetscBool flg;
-    PetscInt  i,nw,nb;
+    PetscInt  i,nw,nl;
     IGABasisType btype[3] = {IGA_BASIS_BSPLINE,IGA_BASIS_BSPLINE,IGA_BASIS_BSPLINE};
-    PetscBool wraps[3]    = {PETSC_FALSE,PETSC_FALSE,PETSC_FALSE};
+    PetscBool    wraps[3] = {PETSC_FALSE, PETSC_FALSE, PETSC_FALSE };
     PetscInt  np,procs[3] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
-    PetscInt  nd,degrs[3] = {2,2,2};
     PetscInt  nq,quadr[3] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
+    PetscInt  ne,elems[3] = {16,16,16};
+    PetscInt  nd,degrs[3] = { 2, 2, 2};
     PetscInt  nc,conts[3] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
-    PetscInt  ne,elems[3] = {16,16,16};
-    PetscReal bbox[3][2]  = {{0,1},{0,1},{0,1}};
+    PetscReal ulims[3][2] = {{0,1},{0,1},{0,1}};
     char      filename[PETSC_MAX_PATH_LEN] = {0};
-    char      vtype[256] = VECSTANDARD;
-    char      mtype[256] = MATBAIJ;
+    char      vtype[256]  = VECSTANDARD;
+    char      mtype[256]  = MATBAIJ;
     PetscInt  dim = (iga->dim > 0) ? iga->dim : 3;
     PetscInt  dof = (iga->dof > 0) ? iga->dof : 1;
     PetscInt  order = iga->order;
 
-    /* Periodicity, degree, and quadrature are initially what they are intially set to */
-    for (i=0; i<dim; i++) wraps[i] = iga->axis[i]->periodic;
-    for (i=0; i<dim; i++) if (iga->axis[i]->p   > 0) degrs[i] = iga->axis[i]->p;
-    for (i=0; i<dim; i++) if (iga->rule[i]->nqp > 0) quadr[i] = iga->rule[i]->nqp;
-    for (i=0; i<dim; i++) btype[i] = iga->basis[i]->type;
+    for (i=0; i<dim; i++) {
+      procs[i] = iga->proc_sizes[i];
+      wraps[i] = iga->axis[i]->periodic;
+      btype[i] = iga->basis[i]->type;
+      if (iga->rule[i]->nqp > 0)
+        quadr[i] = iga->rule[i]->nqp;
+    }
+    for (i=0; i<dim; i++) {
+      IGAAxis axis = iga->axis[i];
+      if (axis->p > 0)
+        degrs[i] = axis->p;
+      else if (i > 0)
+        degrs[i] = degrs[i-1];
+      if (axis->m > 1) {
+        elems[i]    = axis->nel;
+        ulims[i][0] = axis->U[axis->p];
+        ulims[i][1] = axis->U[axis->m-axis->p];
+      } else if (i > 0) {
+        elems[i]    = elems[i-1];
+        ulims[i][0] = ulims[i-1][0];
+        ulims[i][1] = ulims[i-1][1];
+      }
+    }
 
     ierr = PetscObjectOptionsBegin((PetscObject)iga);CHKERRQ(ierr);
 
 
     /* Periodicity */
     ierr = PetscOptionsBoolArray("-iga_periodic","Periodicity","IGAAxisSetPeriodic",wraps,(nw=dim,&nw),&flg);CHKERRQ(ierr);
+    if (flg && nw==0) for (i=0; i<dim; i++) wraps[i] = PETSC_TRUE;
+    if (flg && nw==1) for (i=1; i<dim; i++) wraps[i] = wraps[0]; /* XXX */
     if (flg) for (i=0; i<dim; i++) {
-        PetscBool w = (i<nw) ? wraps[i] : wraps[0];
-        if (nw == 0) w = PETSC_TRUE;
+        PetscBool w = wraps[i];
         ierr = IGAAxisSetPeriodic(iga->axis[i],w);CHKERRQ(ierr);
       }
 
     /* Basis */
     ierr = PetscOptionsEnum("-iga_basis_type","Basis type","IGABasisSetType",IGABasisTypes,(PetscEnum)btype[0],(PetscEnum*)&btype[0],&flg);CHKERRQ(ierr);
-    for (i=0; i<dim; i++) btype[i] = btype[0]; /* XXX */
+    if (flg) for (i=1; i<dim; i++) btype[i] = btype[0]; /* XXX */
     if (flg) for (i=0; i<dim; i++) {
         ierr = IGABasisSetType(iga->basis[i],btype[i]);CHKERRQ(ierr);
       }
     /* 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);
+      ierr = IGAOptionsReject(prefix,"-iga_limits"    );CHKERRQ(ierr);
       ierr = IGARead(iga,filename);CHKERRQ(ierr);
-      ierr = PetscOptionsReject("-iga_elements",  PETSC_NULL);CHKERRQ(ierr);
-      ierr = PetscOptionsReject("-iga_degree",    PETSC_NULL);CHKERRQ(ierr);
-      ierr = PetscOptionsReject("-iga_continuity",PETSC_NULL);CHKERRQ(ierr);
-      ierr = PetscOptionsReject("-iga_limits",    PETSC_NULL);CHKERRQ(ierr);
     } else { /* set axis details */
       ierr = PetscOptionsIntArray ("-iga_elements",  "Elements",  "IGAAxisInitUniform",elems,(ne=dim,&ne),&flg);CHKERRQ(ierr);
       ierr = PetscOptionsIntArray ("-iga_degree",    "Degree",    "IGAAxisSetDegree",  degrs,(nd=dim,&nd),PETSC_NULL);CHKERRQ(ierr);
       ierr = PetscOptionsIntArray ("-iga_continuity","Continuity","IGAAxisInitUniform",conts,(nc=dim,&nc),PETSC_NULL);CHKERRQ(ierr);
-      ierr = PetscOptionsRealArray("-iga_limits",    "Limits",    "IGAAxisInitUniform",&bbox[0][0],(nb=2*dim,&nb),PETSC_NULL);CHKERRQ(ierr);
-      for (i=0; i<dim; i++) {
+      ierr = PetscOptionsRealArray("-iga_limits",    "Limits",    "IGAAxisInitUniform",&ulims[0][0],(nl=2*dim,&nl),PETSC_NULL);CHKERRQ(ierr);
+      for (nl/=2, i=0; i<dim; i++) {
+        IGAAxis axis = iga->axis[i];
+        PetscBool  w = axis->periodic;
+        PetscInt   p = (i<nd) ? degrs[i] : degrs[0];
         PetscInt   N = (i<ne) ? elems[i] : elems[0];
-        PetscInt   p = (i<nd) ? degrs[i] : degrs[0];
+        PetscReal *L = (i<nl) ? ulims[i] : ulims[0];
         PetscInt   C = (i<nc) ? conts[i] : conts[0];
-        PetscReal *U = (i<nb/2) ? &bbox[i][0] : &bbox[0][0];
-        PetscBool  w = iga->axis[i]->periodic;
-        if (p < 1) p = iga->axis[i]->p; if (p < 1) p = 2;
-        if (flg || (iga->axis[i]->p==0||iga->axis[i]->m==1)) {
-          ierr = IGAAxisReset(iga->axis[i]);CHKERRQ(ierr);
-          ierr = IGAAxisSetPeriodic(iga->axis[i],w);CHKERRQ(ierr);
-          ierr = IGAAxisSetDegree(iga->axis[i],p);CHKERRQ(ierr);
-          ierr = IGAAxisInitUniform(iga->axis[i],N,U[0],U[1],C);CHKERRQ(ierr);
+        if (p < 1) {if (axis->p > 0) p = axis->p;   else p =  2;}
+        if (N < 1) {if (axis->m > 1) N = axis->nel; else N = 16;}
+        if (flg || (axis->p==0||axis->m==1)) {
+          ierr = IGAAxisReset(axis);CHKERRQ(ierr);
+          ierr = IGAAxisSetPeriodic(axis,w);CHKERRQ(ierr);
+          ierr = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
+          ierr = IGAAxisInitUniform(axis,N,L[0],L[1],C);CHKERRQ(ierr);
         }
       }
     }
 
-    /* Order */
-    ierr = PetscOptionsInt("-iga_order","Order","IGASetOrder",order,&order,&flg);CHKERRQ(ierr);
-    if (flg) { ierr = IGASetOrder(iga,order);CHKERRQ(ierr);}
-
     /* Quadrature */
+    for (i=0; i<dim; i++) if (quadr[i] < 1) quadr[i] = iga->axis[i]->p + 1;
     ierr = PetscOptionsIntArray("-iga_quadrature","Quadrature points","IGARuleInit",quadr,(nq=dim,&nq),&flg);CHKERRQ(ierr);
     if (flg) for (i=0; i<dim; i++) {
         PetscInt q = (i<nq) ? quadr[i] : quadr[0];
         if (q > 0) {ierr = IGARuleInit(iga->rule[i],q);CHKERRQ(ierr);}
       }
 
+    /* Order */
+    if (order < 0) for (i=0; i<dim; i++) order = PetscMax(order,iga->axis[i]->p);
+    order = PetscMax(order,1); order = PetscMin(order,3);
+    ierr = PetscOptionsInt("-iga_order","Order","IGASetOrder",order,&order,&flg);CHKERRQ(ierr);
+    if (flg) { ierr = IGASetOrder(iga,order);CHKERRQ(ierr);}
+
   setupcalled:
 
     /* Collocation */ {
     if (iga->dof == 1) {ierr = PetscStrcpy(mtype,MATAIJ);CHKERRQ(ierr);}
     if (iga->vectype)  {ierr = PetscStrncpy(vtype,iga->vectype,sizeof(vtype));CHKERRQ(ierr);}
     if (iga->mattype)  {ierr = PetscStrncpy(mtype,iga->mattype,sizeof(mtype));CHKERRQ(ierr);}
-    ierr = PetscOptionsList("-iga_vec_type","Vector type","IGASetVecType",VecList,vtype,vtype,sizeof vtype,&flg);CHKERRQ(ierr);
+    ierr = PetscOptionsList("-iga_vec_type","Vector type","IGASetVecType",VecList,vtype,vtype,sizeof(vtype),&flg);CHKERRQ(ierr);
     if (flg) {ierr = IGASetVecType(iga,vtype);CHKERRQ(ierr);}
-    ierr = PetscOptionsList("-iga_mat_type","Matrix type","IGASetMatType",MatList,mtype,mtype,sizeof mtype,&flg);CHKERRQ(ierr);
+    ierr = PetscOptionsList("-iga_mat_type","Matrix type","IGASetMatType",MatList,mtype,mtype,sizeof(mtype),&flg);CHKERRQ(ierr);
     if (flg) {ierr = IGASetMatType(iga,mtype);CHKERRQ(ierr);}
 
     /* View options, handled in IGASetUp() */
 }
 
 #undef  __FUNCT__
+#define __FUNCT__ "IGAOptionsAlias"
+PetscErrorCode IGAOptionsAlias(const char name[],const char defval[],const char alias[])
+{
+  char           value[4096]= {0};
+  PetscBool      flag = PETSC_FALSE;
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidCharPointer(name,1);
+  PetscValidCharPointer(alias,3);
+  ierr = PetscOptionsHasName(NULL,name,&flag);CHKERRQ(ierr);
+  if (flag) {
+    ierr = PetscOptionsGetString(NULL,name,value,sizeof(value),&flag);CHKERRQ(ierr);
+    ierr = PetscOptionsSetValue(alias,value);CHKERRQ(ierr);
+  } else if (defval) {
+    ierr = PetscOptionsHasName(NULL,alias,&flag);CHKERRQ(ierr);
+    if (!flag) {ierr = PetscOptionsSetValue(alias,defval);CHKERRQ(ierr);}
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
 #define __FUNCT__ "IGACreateSubComms1D"
 PetscErrorCode IGACreateSubComms1D(IGA iga,MPI_Comm subcomms[])
 {
   PetscFunctionReturn(0);
 }
 
+EXTERN_C_BEGIN
+extern PetscReal IGA_Greville(PetscInt i,PetscInt p,const PetscReal U[]);
+extern PetscInt  IGA_FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal U[]);
+EXTERN_C_END
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGASetUp_Stage1"
 PetscErrorCode IGASetUp_Stage1(IGA iga)
 {
-  PetscInt       i;
+  PetscInt       i,dim;
+  PetscInt       grid_sizes[3] = {1,1,1};
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidHeaderSpecific(iga,IGA_CLASSID,1);
   if (iga->setupstage >= 1) PetscFunctionReturn(0);
   iga->setupstage = 1;
 
-  for (i=0; i<iga->dim; i++)
+  dim = iga->dim;
+
+  for (i=0; i<dim; i++)
     {ierr = IGAAxisSetUp(iga->axis[i]);CHKERRQ(ierr);}
-  for (i=iga->dim; i<3; i++)
+  for (i=dim; i<3; i++)
     {ierr = IGAAxisReset(iga->axis[i]);CHKERRQ(ierr);}
 
+  if (!iga->collocation)
+    for (i=0; i<dim; i++)
+      grid_sizes[i] = iga->axis[i]->nel;
+  else
+    for (i=0; i<dim; i++)
+      grid_sizes[i] = iga->axis[i]->nnp;
+
   { /* processor grid and coordinates */
     MPI_Comm    comm = ((PetscObject)iga)->comm;
     PetscMPIInt size,rank;
-    PetscInt    grid_sizes[3] = {1,1,1};
     PetscInt    *proc_sizes = iga->proc_sizes;
     PetscInt    *proc_ranks = iga->proc_ranks;
     ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
     ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
-    for (i=0; i<iga->dim; i++)
-      grid_sizes[i] = iga->axis[i]->nel;
     ierr = IGA_Partition(size,rank,iga->dim,grid_sizes,
                          proc_sizes,proc_ranks);CHKERRQ(ierr);
-    for (i=iga->dim; i<3; i++) {
+    for (i=dim; i<3; i++) {
       proc_sizes[i] = 1;
       proc_ranks[i] = 0;
     }
   }
+
   { /* element partitioning */
     PetscInt *elem_sizes = iga->elem_sizes;
     PetscInt *elem_start = iga->elem_start;
     PetscInt *elem_width = iga->elem_width;
-    for (i=0; i<iga->dim; i++)
-      elem_sizes[i] = iga->axis[i]->nel;
+    for (i=0; i<dim; i++) elem_sizes[i] = grid_sizes[i];
     ierr = IGA_Distribute(iga->dim,iga->proc_sizes,iga->proc_ranks,
                           elem_sizes,elem_width,elem_start);CHKERRQ(ierr);
-    for (i=iga->dim; i<3; i++) {
+    for (i=dim; i<3; i++) {
       elem_sizes[i] = 1;
       elem_start[i] = 0;
       elem_width[i] = 1;
     }
   }
+
+  if (!iga->collocation)
   { /* geometry partitioning */
     PetscInt *geom_sizes  = iga->geom_sizes;
     PetscInt *geom_lstart = iga->geom_lstart;
     PetscInt *geom_lwidth = iga->geom_lwidth;
     PetscInt *geom_gstart = iga->geom_gstart;
     PetscInt *geom_gwidth = iga->geom_gwidth;
-    for (i=0; i<iga->dim; i++) {
+    for (i=0; i<dim; i++) {
       PetscInt nel    = iga->elem_sizes[i];
       PetscInt efirst = iga->elem_start[i];
       PetscInt elast  = iga->elem_start[i] + iga->elem_width[i] - 1;
       geom_gstart[i] = gstart;
       geom_gwidth[i] = gend - gstart;
     }
-    for (i=iga->dim; i<3; i++) {
+    for (i=dim; i<3; i++) {
+      geom_sizes[i]  = 1;
+      geom_lstart[i] = 0;
+      geom_lwidth[i] = 1;
+      geom_gstart[i] = 0;
+      geom_gwidth[i] = 1;
+    }
+  } else
+  {  /* geometry partitioning */
+    PetscInt *geom_sizes  = iga->geom_sizes;
+    PetscInt *geom_lstart = iga->geom_lstart;
+    PetscInt *geom_lwidth = iga->geom_lwidth;
+    PetscInt *geom_gstart = iga->geom_gstart;
+    PetscInt *geom_gwidth = iga->geom_gwidth;
+    for (i=0; i<dim; i++) {
+      PetscInt   p = iga->axis[i]->p;
+      PetscInt   m = iga->axis[i]->m;
+      PetscReal *U = iga->axis[i]->U;
+      PetscInt   n = m - p - 1;
+      geom_sizes[i]  = iga->elem_sizes[i];
+      geom_lstart[i] = iga->elem_start[i];
+      geom_lwidth[i] = iga->elem_width[i];
+      {
+        PetscInt  a = geom_lstart[i];
+        PetscReal u = IGA_Greville(a,p,U);
+        PetscInt  k = IGA_FindSpan(n,p,u,U);
+        geom_gstart[i] = k - p;
+      }
+      {
+        PetscInt  a = geom_lstart[i] + geom_lwidth[i] - 1;
+        PetscReal u = IGA_Greville(a,p,U);
+        PetscInt  k = IGA_FindSpan(n,p,u,U);
+        geom_gwidth[i] = k + 1 - geom_gstart[i];
+      }
+    }
+    for (i=dim; i<3; i++) {
       geom_sizes[i]  = 1;
       geom_lstart[i] = 0;
       geom_lwidth[i] = 1;
       geom_gwidth[i] = 1;
     }
   }
+
   /* element */
   ierr = DMDestroy(&iga->elem_dm);CHKERRQ(ierr);
   /* geometry */
       node_lwidth[i] = iga->geom_lwidth[i];
       node_gstart[i] = iga->geom_gstart[i];
       node_gwidth[i] = iga->geom_gwidth[i];
+      if (rank == 0 && node_lstart[i] < 0) {
+        node_lwidth[i] += node_lstart[i];
+        node_lstart[i] = 0;
+      }
       if (rank == size-1)
         node_lwidth[i] = node_sizes[i] - node_lstart[i];
     }

File src/petiga2d.f90

View file
      if (side==1) k=m(2)
      Xx => Cx(:,:,k); Xw => Cw(:,k)
   end select
-  detJ = 1.0
-  dS = 0.0
+  detJ = 1
+  dS = 0
   do q=1,nqp
      N0(  :) = N(0,:,q)
      N1(1,:) = N(1,:,q)

File src/petiga3d.f90

View file
      if (side==1) k=m(3)
      Xx => Cx(:,:,:,k); Xw => Cw(:,:,k)
   end select
-  detJ = 1.0
-  dS = 0.0
+  detJ = 1
+  dS = 0
   do jq=1,jnqp
      do iq=1,inqp
         do ja=1,jnen; do ia=1,inen

File src/petigaaxis.c

View file
   /* */
   axis->periodic = PETSC_FALSE;
   /* */
+  ierr = PetscMalloc1(2,PetscReal,&axis->U);CHKERRQ(ierr);
   axis->p = 0;
   axis->m = 1;
-  ierr = PetscMalloc1(axis->m+1,PetscReal,&axis->U);CHKERRQ(ierr);
   axis->U[0] = -0.5;
   axis->U[1] = +0.5;
   /* */
+  ierr = PetscMalloc1(1,PetscInt,&axis->span);CHKERRQ(ierr);
   axis->nnp = 1;
   axis->nel = 1;
-  ierr = PetscMalloc1(axis->nel,PetscInt,&axis->span);CHKERRQ(ierr);
   axis->span[0] = 0;
 
   PetscFunctionReturn(0);
   if(base == axis) PetscFunctionReturn(0);
 
   axis->periodic = base->periodic;
+
   axis->p = base->p;
   axis->m = base->m;
   ierr = PetscFree(axis->U);CHKERRQ(ierr);
-  ierr = PetscMalloc((axis->m+1)*sizeof(PetscReal),&axis->U);CHKERRQ(ierr);
+  ierr = PetscMalloc1(axis->m+1,PetscReal,&axis->U);CHKERRQ(ierr);
   ierr = PetscMemcpy(axis->U,base->U,(axis->m+1)*sizeof(PetscReal));CHKERRQ(ierr);
 
-  axis->nel = axis->nnp = 0;
+  axis->nnp = base->nnp;
+  axis->nel = base->nel;
   ierr = PetscFree(axis->span);CHKERRQ(ierr);
+  ierr = PetscMalloc1(axis->nel,PetscInt,&axis->span);CHKERRQ(ierr);
+  ierr = PetscMemcpy(axis->span,base->span,axis->nel*sizeof(PetscInt));CHKERRQ(ierr);
 
   PetscFunctionReturn(0);
 }
   PetscValidPointer(axis,1);
   if (p < 1)
     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
-             "Polynomial degree must be greater than zero, got %D",p);
+             "Polynomial degree must be greater than one, got %D",p);
+  if (axis->p > 0 && axis->m > 1 && axis->p != p)
+    SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,
+             "Cannot change degree to %D after it was set to %D",p,axis->p);
   axis->p = p;
   PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAAxisSetKnots"
-PetscErrorCode IGAAxisSetKnots(IGAAxis axis,PetscInt m,PetscReal U[])
+PetscErrorCode IGAAxisSetKnots(IGAAxis axis,PetscInt m,const PetscReal U[])
 {
-  PetscInt       k;
-  PetscReal      *V = 0;
+  PetscInt       p,n,k;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(axis,1);
-  if (U) PetscValidPointer(U,3);
+  PetscValidPointer(U,3);
 
   if (axis->p < 1)
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,
             "Must call IGAAxisSetDegree() first");
   if (m < 2*axis->p+1)
     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
-             "Number of knots must be at least %D, got %D",2*(axis->p+1),m+1);
-  if (U) for (k=1; k<=m; k++)
-           if (U[k-1] > U[k])
-             SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
-                      "Knot sequence must be increasing, "
-                      "got U[%D]=%G > U[%D]=%G",
-                      k-1,U[k-1],k,U[k]);
+             "Number of knots must be at least %D, got %D",
+             2*(axis->p+1),m+1);
 
-  if (m != axis->m || !axis->U) {
-    ierr = PetscMalloc((m+1)*sizeof(PetscReal),&V);CHKERRQ(ierr);
-    ierr = PetscMemzero(V,(m+1)*sizeof(PetscReal));CHKERRQ(ierr);
+  p = axis->p;
+  n = m -p - 1;
+
+  for (k=1; k<=m;) {
+    PetscInt i = k, s = 1;
+    /* check increasing sequence */
+    if (U[k-1] > U[k])
+      SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
+               "Knot sequence must be increasing, got U[%D]=%G > U[%D]=%G",
+               k-1,U[k-1],k,U[k]);
+    /* check multiplicity */
+    while (++k < m && U[i] == U[k]) s++;
+    if (s > p)
+      SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
+               "Knot U[%D]=%G has multiplicity %D greater than polynomial degree %D",
+               i,U[i],s,p);
+  }
+
+  if (m != axis->m) {
+    PetscReal *V;
+    ierr = PetscMalloc1(m+1,PetscReal,&V);CHKERRQ(ierr);
     ierr = PetscFree(axis->U);CHKERRQ(ierr);
     axis->m = m;
     axis->U = V;
   }
-  if (U) { ierr = PetscMemcpy(axis->U,U,(m+1)*sizeof(PetscReal));CHKERRQ(ierr); }
+  ierr = PetscMemcpy(axis->U,U,(m+1)*sizeof(PetscReal));CHKERRQ(ierr);
 
-  axis->nel = axis->nnp = 0;
+  axis->nel = 0;
   ierr = PetscFree(axis->span);CHKERRQ(ierr);
+  ierr = IGAAxisGetSpans(axis,&axis->nel,&axis->span);CHKERRQ(ierr);
+
+  if (axis->periodic) {
+    PetscInt s = 1;
+    while(s < p && U[m-p] == U[m-p+s]) s++;
+    axis->nnp = n-p+s;
+  } else {
+    axis->nnp = n+1;
+  }
 
   PetscFunctionReturn(0);
 }
   PetscFunctionReturn(0);
 }
 
+EXTERN_C_BEGIN
+extern PetscInt IGA_SpanCount(PetscInt n,PetscInt p,const PetscReal U[]);
+extern PetscInt IGA_SpanIndex(PetscInt n,PetscInt p,const PetscReal U[],PetscInt index[]);
+EXTERN_C_END
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAAxisGetSpans"
+PetscErrorCode IGAAxisGetSpans(IGAAxis axis,PetscInt *nel,PetscInt *span[])
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(axis,1);
+  if (nel)  PetscValidIntPointer(nel,2);
+  if (span) PetscValidPointer(span,3);
+  if (!axis->span) {
+    PetscInt p = axis->p;
+    PetscInt m = axis->m;
+    PetscInt n = m - p - 1;
+    axis->nel = IGA_SpanCount(n,p,axis->U);
+    ierr = PetscMalloc1(axis->nel,PetscInt,&axis->span);CHKERRQ(ierr);
+    (void)IGA_SpanIndex(n,p,axis->U,axis->span);
+  }
+  if (nel)  *nel  = axis->nel;
+  if (span) *span = axis->span;
+  PetscFunctionReturn(0);
+}
+
+#undef  __FUNCT__
+#define __FUNCT__ "IGAAxisInit"
+PetscErrorCode IGAAxisInit(IGAAxis axis,PetscInt p,PetscInt m,const PetscReal U[])
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidPointer(axis,1);
+  axis->p = 0;
+  ierr = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
+  ierr = IGAAxisSetKnots(axis,m,U);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
 #undef  __FUNCT__
 #define __FUNCT__ "IGAAxisInitBreaks"
 PetscErrorCode IGAAxisInitBreaks(IGAAxis axis,PetscInt nu,const PetscReal u[],PetscInt C)
     ierr = PetscFree(axis->U);CHKERRQ(ierr);
     axis->m = m;
     axis->U = U;
-  } else U = axis->U;
+  } else
+    U = axis->U;
 
   for(k=0; k<=p; k++) { /* open part */
     U[k]   = u[0];
     }
   }
 
-  axis->nel = axis->nnp = 0;
+  axis->nel = r;
   ierr = PetscFree(axis->span);CHKERRQ(ierr);
+  ierr = PetscMalloc1(axis->nel,PetscInt,&axis->span);CHKERRQ(ierr);
+  for(i=0; i<axis->nel; i++) axis->span[i] = p + i*s;
+
+  axis->nnp = axis->periodic ? n-C : n+1;
 
   PetscFunctionReturn(0);
 }
     ierr = PetscFree(axis->U);CHKERRQ(ierr);
     axis->m = m;
     axis->U = U;
-  } else U = axis->U;
+  } else
+    U = axis->U;
 
   for(k=0; k<=p; k++) { /* open part */
     U[k]   = Ui;
     }
   }
 
-  axis->nel = axis->nnp = 0;
+  axis->nel = r;
   ierr = PetscFree(axis->span);CHKERRQ(ierr);
+  ierr = PetscMalloc1(axis->nel,PetscInt,&axis->span);CHKERRQ(ierr);
+  for(i=0; i<axis->nel; i++) axis->span[i] = p + i*s;
+
+  axis->nnp = axis->periodic ? n-C : n+1;
 
   PetscFunctionReturn(0);
 }
 
-EXTERN_C_BEGIN
-extern PetscInt IGA_SpanCount(PetscInt n,PetscInt p,const PetscReal U[]);
-extern PetscInt IGA_SpanIndex(PetscInt n,PetscInt p,const PetscReal U[],PetscInt index[]);
-EXTERN_C_END
-
 #undef  __FUNCT__
 #define __FUNCT__ "IGAAxisSetUp"
 PetscErrorCode IGAAxisSetUp(IGAAxis axis)
   if (axis->p < 1)
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,
             "Must call IGAAxisSetDegree() first");
-  if (!axis->U || axis->m < 2*axis->p+1)
+  if (axis->m < 2*axis->p+1)
     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,
             "Must call IGAAxisSetKnots() first");
 
   n = m - p - 1;
   U = axis->U;
 
-#ifdef PETSC_USE_DEBUG
-  {
-    PetscInt k = 1;
-    while (k <= m) {
-      PetscInt i = k, s = 1;
-      /* check increasing sequence */
-      if (U[k-1] > U[k])
-        SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
-                 "Knot sequence must be increasing, "
-                 "got U[%D]=%G > U[%D]=%G",
-                 k-1,U[k-1],k,U[k]);
-      /* check multiplicity */
-      while (++k < m && U[k-1] == U[k]) s++;
-      if (s > p)
-        SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,
-                 "Knot U[%D]=%G has multiplicity %D "
-                 "greater than polynomial degree %D",
-                 i,U[i],s,p);
-    }
-  }
-#endif
-
-  axis->nel = axis->nnp = 0;
-  ierr = PetscFree(axis->span);CHKERRQ(ierr);
-
-  axis->nel = IGA_SpanCount(n,p,U);
-  ierr = PetscMalloc1(axis->nel,PetscInt,&axis->span);CHKERRQ(ierr);
-  (void)IGA_SpanIndex(n,p,U,axis->span);
+  ierr = IGAAxisGetSpans(axis,&axis->nel,&axis->span);CHKERRQ(ierr);
 
   if (axis->periodic) {
     PetscInt s = 1;

File src/petigabasis.c

View file
                       "Derivative order must be grather than zero, got %D",d);
 
   p = axis->p;
-  n = axis->m - p -1;
+  n = axis->m - p - 1;
   U = axis->U;
 
   nel  = axis->nnp;
     PetscInt  k = IGA_FindSpan(n,p,u,U);
     PetscReal *N = &value[iel*nen*ndr];
     offset[iel] = k-p;
+    detJ[iel]   = U[k+1]-U[k];
     point[iel]  = u;
-    detJ[iel]   = U[k+1]-U[k];
     IGA_Basis_BSpline(k,u,p,d,U,N);
   }
 
   basis->point  = point;
   basis->value  = value;
 
+  {
+    PetscInt  k0 = p,    k1 = n;
+    PetscReal u0 = U[p], u1 = U[n+1];
+    ierr = PetscMalloc1(nen*ndr,PetscReal,&basis->bnd_value[0]);CHKERRQ(ierr);
+    ierr = PetscMalloc1(nen*ndr,PetscReal,&basis->bnd_value[1]);CHKERRQ(ierr);
+    basis->bnd_offset[0] = k0-p; basis->bnd_offset[1] =  k1-p;
+    basis->bnd_detJ  [0] =  1.0; basis->bnd_detJ  [1] =   1.0;
+    basis->bnd_weight[0] =  1.0; basis->bnd_weight[1] =   1.0;
+    basis->bnd_point [0] =   u0; basis->bnd_point [1] =    u1;
+    IGA_Basis_BSpline(k0,u0,p,d,U,basis->bnd_value[0]);
+    IGA_Basis_BSpline(k1,u1,p,d,U,basis->bnd_value[1]);
+  }
+
   PetscFunctionReturn(0);
 }
 
 PetscInt IGA_FindSpan(PetscInt n,PetscInt p,PetscReal u, const PetscReal U[])
 {
   PetscInt low,high,span;
+  if(u <= U[p])   return p;
   if(u >= U[n+1]) return n;
-  if(u <= U[p])   return p;
   low  = p;
   high = n+1;
   span = (high+low)/2;

File src/petigabsp.f90

View file
   real   (kind=IGA_REAL_KIND   )  :: saved, temp, d
   real   (kind=IGA_REAL_KIND   )  :: left(p), right(p)
   real   (kind=IGA_REAL_KIND   )  :: ndu(0:p,0:p), a(0:1,0:p)
-  ndu(0,0) = 1.0
+  ndu(0,0) = 1
   do j = 1, p
      left(j)  = uu - U(i+1-j)
      right(j) = U(i+j) - uu
-     saved = 0.0
+     saved = 0
      do r = 0, j-1
         ndu(j,r) = right(r+1) + left(j-r)
         temp = ndu(r,j-1) / ndu(j,r)
   ders(:,0) = ndu(:,p)
   do r = 0, p
      s1 = 0; s2 = 1;
-     a(0,0) = 1.0
+     a(0,0) = 1
      do k = 1, n
-        d = 0.0
+        d = 0
         rk = r-k; pk = p-k;
         if (r >= k) then
            a(s2,0) = a(s1,0) / ndu(pk+1,rk)
   end do
 
   do m = 0, p
-     Lp = 1.0
+     Lp = 1
      do i = 0, p
         if (i == m) cycle
         Lp = Lp * (uu-X(i))/(X(m)-X(i))
 
   if (d < 1) return
   do m = 0, p
-     Ls1 = 0.0
+     Ls1 = 0
      do j = 0, p
         if (j == m) cycle
-        Lp = 1.0
+        Lp = 1
         do i = 0, p
            if (i == m .or. i == j) cycle
            Lp = Lp * (uu-X(i))/(X(m)-X(i))
 
   if (d < 2) return
   do m = 0, p
-     Ls2 = 0.0
+     Ls2 = 0
      do k = 0, p
         if (k == m) cycle
-        Ls1 = 0.0
+        Ls1 = 0
         do j = 0, p
            if (j == m .or. j == k) cycle
-           Lp = 1.0
+           Lp = 1
            do i = 0, p
               if (i == m .or. i == k .or. i == j) cycle
               Lp = Lp * (uu-X(i))/(X(m)-X(i))
 
   if (d < 3) return
   do m = 0, p
-     Ls3 = 0.0
+     Ls3 = 0
      do l = 0, p
         if (l == m) cycle
-        Ls2 = 0.0
+        Ls2 = 0
         do k = 0, p
            if (k == m .or. k == l) cycle
-           Ls1 = 0.0
+           Ls1 = 0
            do j = 0, p
               if (j == m .or. j == l .or. j == k) cycle
-              Lp = 1.0
+              Lp = 1
               do i = 0, p
                  if (i == m .or. i == l .or. i == k .or. i == j) cycle
                  Lp = Lp * (uu-X(i))/(X(m)-X(i))
   real   (kind=IGA_REAL_KIND   ), intent(out)      :: B(0:d,0:p)
   integer(kind=IGA_INTEGER_KIND)  :: i, k
   real   (kind=IGA_REAL_KIND   )  :: J, x, Lp(0:p,0:d)
-  real   (kind=IGA_REAL_KIND   ), parameter :: two = 2.0
+  real   (kind=IGA_REAL_KIND   ), parameter :: two = 2
 
-  J = (U(kk+1)-U(kk))/2.0
-  x = (uu-U(kk))/J - 1.0
+  J = (U(kk+1)-U(kk))/2
+  x = (uu-U(kk))/J - 1
 
-  B(0,0) = (1.0-x)/2.0
-  B(0,p) = (x+1.0)/2.0
+  B(0,0) = (1-x)/2
+  B(0,p) = (x+1)/2
   if (d > 0) then
-     B(1,0) = -0.5
-     B(1,p) = +0.5
+     B(1,0) = -1/two
+     B(1,p) = +1/two
   endif
 
   if (p > 1) then
-     Lp(:,:) = 0.0
-     Lp(0,0) = 1.0
+     Lp(:,:) = 0
+     Lp(0,0) = 1
      Lp(1,0) = x
      if (d > 0) then
-        Lp(0,1) = 0.0
-        Lp(1,1) = 1.0
+        Lp(0,1) = 0
+        Lp(1,1) = 1
      end if
      do i = 1, p-1
         Lp(i+1,0) = ((2*i+1)*x*Lp(i,0) - i*Lp(i-1,0))/(i+1)

File src/petigaelem.c

View file
   element->count =  0;
   element->index = -1;
 
+  if (element->rowmap != element->mapping)
+    {ierr = PetscFree(element->rowmap);CHKERRQ(ierr);}
+  if (element->colmap != element->mapping)
+    {ierr = PetscFree(element->colmap);CHKERRQ(ierr);}
+  element->rowmap = element->colmap = 0;
   ierr = PetscFree(element->mapping);CHKERRQ(ierr);
+
   ierr = PetscFree(element->geometryX);CHKERRQ(ierr);
   ierr = PetscFree(element->rationalW);CHKERRQ(ierr);
   ierr = PetscFree(element->propertyA);CHKERRQ(ierr);
   PetscInt *start;
   PetscInt *width;
   PetscInt *sizes;
+  IGABasis *BD;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(element,1);
   start = iga->elem_start;
   width = iga->elem_width;
   sizes = iga->elem_sizes;
-  element->BD = iga->basis;
-  if (PetscUnlikely(element->collocation)) { /* collocation */
-    start = iga->node_lstart;
-    width = iga->node_lwidth;
-    sizes = iga->node_sizes;
-    element->BD = iga->node_basis;
-  }
+  BD = element->collocation ? iga->node_basis : iga->basis;
+
   { /* */
     PetscInt i,dim = element->dim;
     PetscInt nel=1,nen=1,nqp=1;
+    for (i=0; i<3; i++)
+      element->BD[i] = BD[i];
     for (i=0; i<dim; i++) {
       element->start[i] = start[i];
       element->width[i] = width[i];
       element->sizes[i] = sizes[i];
-      nel *= element->width[i];
-      nen *= element->BD[i]->nen;
-      nqp *= element->BD[i]->nqp;
+      nel *= width[i];
+      nen *= BD[i]->nen;
+      nqp *= BD[i]->nqp;
     }
     for (i=dim; i<3; i++) {
       element->start[i] = 0;
     element->nen   = nen;
     element->nqp   = nqp;
   }
+  { /**/
+    PetscInt nen = element->nen;
+    ierr = PetscMalloc1(nen,PetscInt,&element->mapping);CHKERRQ(ierr);
+    if (!element->collocation) {
+      element->neq = nen;
+      element->rowmap = element->mapping;
+    } else {
+      element->neq = 1;
+      ierr = PetscMalloc1(1,PetscInt,&element->rowmap);CHKERRQ(ierr);
+    }
+    element->colmap = element->mapping;
+  }
   { /* */
     PetscInt dim = element->dim;
     PetscInt nsd = element->nsd;
 
     /* */
 
-    ierr = PetscMalloc1(nen,PetscInt,&element->mapping);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*nsd,PetscReal,&element->geometryX);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen,PetscReal,&element->rationalW);CHKERRQ(ierr);
     ierr = PetscMalloc1(nen*npd,PetscScalar,&element->propertyA);CHKERRQ(ierr);
 
     ierr = PetscMemzero(element->normal,sizeof(PetscReal)*nqp*dim);CHKERRQ(ierr);
   }
-  element->neq    = element->nen;
-  element->rowmap = element->mapping;
-  element->colmap = element->mapping;
-  if (PetscUnlikely(element->collocation)) { /* collocation */
-    element->neq    = 1;
-    element->rowmap = &element->index;
-    element->colmap = element->mapping;
-  }
   { /* */
     PetscInt nen = element->nen;
     PetscInt dof = element->dof;
   element->nval = 0;
   element->nvec = 0;
   element->nmat = 0;
+  element->boundary_id = -1;
   /* */
   index = ++element->index;
   if (PetscUnlikely(index >= element->count)) {
 #define __FUNCT__ "IGAElementNextUserOps"
 PetscBool IGAElementNextUserOps(IGAElement element,IGAUserOps *userops)
 {
-  static PetscInt counter = 0; /* XXX */
-  IGA iga = element->parent;
+  IGA      iga = element->parent;
   PetscInt dim = element->dim;
-  for (; PetscUnlikely(counter < 2*dim); counter++) {
-    PetscInt dir  = counter / 2;
-    PetscInt side = counter % 2;
-    PetscInt iel  = side ? element->sizes[dir]-1 : 0;
-    if (element->ID[dir] != iel) continue;
+  while (++element->boundary_id < 2*dim) {
+    PetscInt i = element->boundary_id / 2;
+    PetscInt s = element->boundary_id % 2;
+    PetscInt e = s ? element->sizes[i]-1 : 0;
+    if (element->ID[i] != e) continue;
+    if (!iga->boundary[i][s]->userops) continue;
+    *userops = iga->boundary[i][s]->userops;
     element->atboundary = PETSC_TRUE;
-    element->atboundary_id = counter;
-    *userops = iga->boundary[dir][side]->userops;
-    if (!(*userops)) continue;
-    counter++;
     return PETSC_TRUE;
   }
-  if (PetscLikely(counter == 2*dim)) {
+  if (element->boundary_id++ == 2*dim) {
+    *userops = iga->userops;
     element->atboundary = PETSC_FALSE;
-    element->atboundary_id = 0;
-    *userops = iga->userops;
-    counter++;
     return PETSC_TRUE;
   }
   *userops = 0;
-  counter = 0;
+  element->atboundary  = PETSC_FALSE;
+  element->boundary_id = -1;
   return PETSC_FALSE;
 }
 
   (*point)->nsd = element->nsd;
   (*point)->npd = element->npd;
 
-  if (PetscUnlikely(!element->atboundary)) {
+  if (PetscLikely(!element->atboundary)) {
     ierr = IGAElementBuildQuadrature(element);CHKERRQ(ierr);
     ierr = IGAElementBuildShapeFuns(element);CHKERRQ(ierr);
   } else {
-    PetscInt i = element->atboundary_id / 2;
-    PetscInt s = element->atboundary_id % 2;
+    PetscInt i = element->boundary_id / 2;
+    PetscInt s = element->boundary_id % 2;
     (*point)->count = element->nqp / element->BD[i]->nqp;
     ierr = IGAElementBuildQuadratureAtBoundary(element,i,s);CHKERRQ(ierr);
     ierr = IGAElementBuildShapeFunsAtBoundary (element,i,s);CHKERRQ(ierr);
-    PetscFunctionReturn(0);
   }
   PetscFunctionReturn(0);
 }
     PetscFunctionReturn(PETSC_ERR_PLIB);
   }
   *point = PETSC_NULL;
-  PetscFunctionReturn(0);
+  /* XXX */
+  if (PetscLikely(!element->collocation)) PetscFunctionReturn(0);
+  if (PetscLikely(!element->atboundary))  PetscFunctionReturn(0);
+  element->atboundary  = PETSC_FALSE;
+  element->boundary_id = 2*element->dim;
+  /* XXX */
+ PetscFunctionReturn(0);
 }
 
 #undef  __FUNCT__
         }
       }
     }
+    if (PetscUnlikely(element->collocation)) {
+      PetscInt iA = ID[0] - istart;
+      PetscInt jA = ID[1] - jstart;
+      PetscInt kA = ID[2] - kstart;
+      element->rowmap[0] = iA + jA*jstride + kA*kstride;
+    }
   }
   PetscFunctionReturn(0);
 }
   }
 }
 
+PETSC_STATIC_INLINE
+IGABoundary AtBoundary(IGAElement element,PetscInt dir,PetscInt side)
+{
+  IGABoundary b = element->parent->boundary[dir][side];
+  PetscInt e = side ? element->sizes[dir]-1 : 0;
+  return (element->ID[dir] == e) ? b : PETSC_NULL;
+}
+
+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__ "IGAElementBuildFix"
 PetscErrorCode IGAElementBuildFix(IGAElement element)
       if (ID[i] == 0 && !w) BuildFix(element,i,0);
       if (ID[i] == e && !w) BuildFix(element,i,1);
     }
-  } 
+  }
   PetscFunctionReturn(0);
- collocation: 
+ collocation:
   element->nfix  = 0;
   element->nflux = 0;
   {
-    PetscInt A0[3] = {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT};
-    PetscInt A1[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
+    PetscInt L[3] = {PETSC_MIN_INT,PETSC_MIN_INT,PETSC_MIN_INT};
+    PetscInt R[3] = {PETSC_MAX_INT,PETSC_MAX_INT,PETSC_MAX_INT};
     {
       IGAAxis  *AX = element->parent->axis;
       PetscInt i,dim = element->dim;
       for (i=0; i<dim; i++) {
         PetscBool w = AX[i]->periodic;
         PetscInt  n = element->sizes[i]-1; /* last node */
-        A0[i] = 0; if (!w) A1[i] = n;
+        L[i] = 0; if (!w) R[i] = n;
       }
     }
     {
               PetscInt iA = ioffset + ia;
               PetscInt jA = joffset + ja;
               PetscInt kA = koffset + ka;
-              /**/ if (iA == A0[0]) AddFixa(element,b[0][0],a);
-              else if (iA == A1[0]) AddFixa(element,b[0][1],a);
-              /**/ if (jA == A0[1]) AddFixa(element,b[1][0],a);
-              else if (jA == A1[1]) AddFixa(element,b[1][1],a);
-              /**/ if (kA == A0[2]) AddFixa(element,b[2][0],a);
-              else if (kA == A1[2]) AddFixa(element,b[2][1],a);
+              /**/ if (iA == L[0]) AddFixa(element,b[0][0],a);
+              else if (iA == R[0]) AddFixa(element,b[0][1],a);
+              /**/ if (jA == L[1]) AddFixa(element,b[1][0],a);
+              else if (jA == R[1]) AddFixa(element,b[1][1],a);
+              /**/ if (kA == L[2]) AddFixa(element,b[2][0],a);
+              else if (kA == R[2]) AddFixa(element,b[2][1],a);
               a++;
             }
     }
   PetscFunctionReturn(0);
  collocation:
   {
-    PetscInt M = element->neq * element->dof;
-    PetscInt N = element->nen * element->dof;
-    PetscInt f,n;
-    n = element->nfix;
-    for (f=0; f<n; f++) {
-      PetscInt k = element->ifix[f];
-      PetscInt a = k / element->dof;
-      if (element->index == element->colmap[a]) {
-        PetscInt    c = k % element->dof;
-        PetscScalar v = element->vfix[f];
-        PetscInt i,j;
-        for (i=0; i<M; i++) K[i*N+k] = 0.0;
-        for (j=0; j<N; j++) K[c*N+j] = 0.0;
-        K[c*N+k] = 1.0;
-        F[c] = v;
+    PetscInt dim = element->dim;
+    PetscInt dof = element->dof;
+    PetscInt nen = element->nen;
+    PetscInt N = nen * dof;
+    PetscInt dir,side;
+    for (dir=0; dir<dim; dir++) {
+      for (side=0; side<2; side++) {
+        IGABoundary b = AtBoundary(element,dir,side);
+        if(b && b->nload) {
+          PetscInt  f,n = b->nload;
+          PetscReal normal[3] = {0.0,0.0,0.0};
+          PetscReal *dshape;
+          if (!element->geometry) {
+            normal[dir] = side ? 1.0 : -1.0;
+            dshape = element->basis[1];
+          } else {
+            PetscReal dS, *F = element->gradX[0];
+            IGA_GetNormal(dim,dir,side,F,&dS,normal);
+            dshape = element->shape[1];
+          }
+          for (f=0; f<n; f++) {
+            PetscInt  c = b->iload[f];
+            PetscReal v = b->vload[f];
+            PetscInt  a,j;
+            for (j=0; j<N; j++) K[c*N+j] = 0.0;
+            for (a=0; a<nen; a++)
+              K[c*N+a*dof+c] = DOT(dim,&dshape[a*dim],normal);
+            F[c] = v;
+          }
+        }
+        if (b && b->count) {
+          PetscInt  f,n = b->count;
+          PetscReal *shape;
+          if (!element->geometry) {
+            shape = element->basis[0];
+          } else {
+            shape = element->shape[0];
+          }
+          for (f=0; f<n; f++) {
+            PetscInt  c = b->field[f];
+            PetscReal v = b->value[f];
+            PetscInt  a,j;
+            for (j=0; j<N; j++) K[c*N+j] = 0.0;
+            for (a=0; a<nen; a++)
+              K[c*N+a*dof+c] = shape[a];
+            F[c] = v;
+          }
+        }
       }
     }
   }
     }
   }
   PetscFunctionReturn(0);
- collocation: 
+ collocation:
   {
     PetscInt f,n;
     n = element->nfix;
     for (f=0; f<n; f++) {
       PetscInt k = element->ifix[f];
       PetscInt a = k / element->dof;
-      if (element->index == element->colmap[a]) {
+      if (element->rowmap[0] == element->colmap[a]) {
         PetscInt    c = k % element->dof;
         PetscScalar v = element->vfix[f];
         PetscScalar u = element->ufix[f];
     for (f=0; f<n; f++) {
       PetscInt k = element->ifix[f];
       PetscInt a = k / element->dof;
-      if (element->index == element->colmap[a]) {
+      if (element->rowmap[0] == element->colmap[a]) {
         PetscInt c = k % element->dof;
         PetscInt i,j;
-        if (element->index != element->colmap[a]) continue;
         for (i=0; i<M; i++) J[i*N+k] = 0.0;
         for (j=0; j<N; j++) J[c*N+j] = 0.0;
         J[c*N+k] = 1.0;

File src/petigaftn.F90

View file
       scalar (kind=IGA_SCALAR_KIND)  :: V
       ! V = dot_product(N,U)
       integer a
-      V = 0.0
+      V = 0
       do a = 1, size(U,1) ! nen
          V = V + N(a) * U(a)
       end do
       scalar (kind=IGA_SCALAR_KIND)  :: V(DOF)            ! dof
       ! V = matmul(N,transpose(U))
       integer a
-      V = 0.0
+      V = 0
       do a = 1, size(U,2) ! nen
          V = V + N(a) * U(:,a)
       end do
       scalar (kind=IGA_SCALAR_KIND)  :: V(DIM)            ! dim
       !V = matmul(N,U)
       integer a
-      V = 0.0
+      V = 0
       do a = 1, size(U,1) ! nen
          V(:) = V(:) + N(:,a) * U(a)
       end do
       scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DOF)        ! dim,dof
       ! V = matmul(N,transpose(U))
       integer a, c
-      V = 0.0
+      V = 0
       do a = 1, size(U,2) ! nen
          do c = 1, size(U,1) ! dof
             V(:,c) = V(:,c) + N(:,a) * U(c,a)
       scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:)     ! nen
       scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM)          ! dim,dim
       integer a
-      V = 0.0
+      V = 0
       do a = 1, size(U,1) ! nen
          V(:,:) = V(:,:) + N(:,:,a) * U(a)
       end do
       scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:)   ! dof,nen
       scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DOF)      ! dim,dim,dof
       integer a, c
-      V = 0.0
+      V = 0
       do a = 1, size(U,2) ! nen
          do c = 1, size(U,1) ! dof
             V(:,:,c) = V(:,:,c) + N(:,:,a) * U(c,a)
       scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:)       ! nen
       scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DIM)        ! dim,dim,dim
       integer a
-      V = 0.0
+      V = 0
       do a = 1, size(U,1) ! nen
          V(:,:,:) = V(:,:,:) + N(:,:,:,a) * U(a)
       end do
       scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:)     ! dof,nen
       scalar (kind=IGA_SCALAR_KIND)  :: V(DIM,DIM,DIM,DOF)    ! dim,dim,dim,dof
       integer a, c
-      V = 0.0
+      V = 0
       do a = 1, size(U,2) ! nen
          do c = 1, size(U,1) ! dof
             V(:,:,:,c) = V(:,:,:,c) + N(:,:,:,a) * U(c,a)
       scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:) ! dim,nen
       scalar (kind=IGA_SCALAR_KIND)  :: V
       integer a, c, i
-      V = 0.0
+      V = 0
       do a = 1, size(U,2) ! nen
          do i = 1, size(N,1) ! dim
             V = V + N(i,a) * U(i,a)
       scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:)     ! nen
       scalar (kind=IGA_SCALAR_KIND)  :: V
       integer a, c, i
-      V = 0.0
+      V = 0
       do a = 1, size(U,1) ! nen
          do i = 1, size(N,1) ! dim
             V = V + N(i,i,a) * U(a)
       scalar (kind=IGA_SCALAR_KIND), intent(in) :: U(:,:)   ! dof,nen
       scalar (kind=IGA_SCALAR_KIND)  :: V(DOF)
       integer a, c, i
-      V = 0.0
+      V = 0
       do a = 1, size(U,2) ! nen
          do c = 1, size(U,1) ! dof
             do i = 1, size(N,1) ! dim

File src/petigageo.f90.in

View file
   ! 1st derivatives
   X1 = transpose(Grad)
   E1 = transpose(InvG)
-  R1 = 0.0
+  R1 = 0
   do idx = 1,nen
      do i = 1,dim
         do a = 1,dim
 
   ! 2nd derivatives
   if (order < 2) return
-  X2 = 0.0
+  X2 = 0
   do b = 1,dim
      do a = 1,dim
         do i = 1,dim
         end do
      end do
   end do
-  E2 = 0.0
+  E2 = 0
   do j = 1,dim
      do i = 1,dim
         do c = 1,dim
         end do
      end do
   end do
-  R2 = 0.0
+  R2 = 0
   do idx = 1,nen
      do j = 1,dim
         do i = 1,dim
 
   ! 3rd derivatives
   if (order < 3) return
-  X3 = 0.0
+  X3 = 0
   do c = 1,dim
      do b = 1,dim
         do a = 1,dim
         end do
      end do
   end do
-  E3 = 0.0
+  E3 = 0
   do k = 1,dim
      do j = 1,dim
         do i = 1,dim
         end do
      end do
   end do
-  R3 = 0.0
+  R3 = 0
   do idx = 1,nen
      do k = 1,dim
         do j = 1,dim

File src/petigainv.f90.in

View file
             - A(2,1) * ( A(1,2)*A(3,3) - A(3,2)*A(1,3) ) &
             + A(3,1) * ( A(1,2)*A(2,3) - A(2,2)*A(1,3) )
   case default
-     detA = 0.0
+     detA = 0
   end select
 end function Determinant
 
   real   (kind=IGA_REAL_KIND   )  :: invA(dim,dim)
   select case (dim)
   case (1)
-     invA = 1.0/detA
+     invA = 1/detA
   case (2)
      invA(1,1) = + A(2,2)
      invA(2,1) = - A(2,1)
      invA(3,3) = + A(1,1)*A(2,2) - A(1,2)*A(2,1)
      invA = invA/detA
   case default
-     invA = 0.0
+     invA = 0
   end select
 end function Inverse

File src/petigaio.c