Commits

Lisandro Dalcin committed e1d8929

Change IGAAxisInitUniform() to required previos call to IGAAxisSetOrder()

  • Participants
  • Parent commits f009588

Comments (0)

Files changed (14)

demo/CahnHilliard2D.c

   IGAAxis axis0;
   ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
   ierr = IGAAxisSetPeriodic(axis0,PETSC_TRUE);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,p,C,N,0.0,1.0);CHKERRQ(ierr);
+  ierr = IGAAxisSetOrder(axis0,p);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,N,0.0,1.0,C);CHKERRQ(ierr);
   IGAAxis axis1;
   ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
   ierr = IGAAxisCopy(axis0,axis1);CHKERRQ(ierr);

demo/InputOutput.c

   if (n3<3) C[2] = C[0]; if (n3<2) C[1] = C[0];
   for (i=0; i<dim; i++)  if (C[i] ==-1) C[i] = p[i] - 1;
   
-  IGA iga; IGAAxis axis;
+  IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
   ierr = IGASetDof(iga,dof);CHKERRQ(ierr);
   for (i=0; i<dim; i++) {
+    IGAAxis axis;
     ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
     ierr = IGAAxisSetPeriodic(axis,b[i]);CHKERRQ(ierr);
-    ierr = IGAAxisInitUniform(axis,p[i],C[i],N[i],0.0,1.0);CHKERRQ(ierr);
+    ierr = IGAAxisSetOrder(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],0.0,1.0,C[i]);CHKERRQ(ierr);
   }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,1);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-  IGAAxis axis0;
-  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,p,C,N,-1.0,1.0);CHKERRQ(ierr);
+  IGAAxis axis;
+  ierr = IGAGetAxis(iga,0,&axis);CHKERRQ(ierr);
+  ierr = IGAAxisSetOrder(axis,p);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis,N,-1.0,1.0,C);CHKERRQ(ierr);
 
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,2);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
-  IGAAxis axis0;
-  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,p[0],C[0],N[0],-1.0,1.0);CHKERRQ(ierr);
-  IGAAxis axis1;
-  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis1,p[1],C[1],N[1],-1.0,1.0);CHKERRQ(ierr);
-
+  PetscInt i;
+  for (i=0; i<2; i++) {
+    IGAAxis axis;
+    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
+  ierr = IGAAxisSetOrder(axis,p[i]);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis,N[i],-1.0,1.0,C[i]);CHKERRQ(ierr);
+  }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
 
   if (n3<3) C[2] = C[0]; if (n3<2) C[1] = C[0];
   for (i=0; i<dim; i++)  if (C[i] ==-1) C[i] = p[i] - 1;
   
-  IGA iga; IGAAxis axis;
+  IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
   ierr = IGASetDof(iga,dof);CHKERRQ(ierr);
   for (i=0; i<dim; i++) {
+    IGAAxis axis;
     ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
     ierr = IGAAxisSetPeriodic(axis,t[i]);CHKERRQ(ierr);
-    ierr = IGAAxisInitUniform(axis,p[i],C[i],N[i],0.0,1.0);CHKERRQ(ierr);
+    ierr = IGAAxisSetOrder(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],0.0,1.0,C[i]);CHKERRQ(ierr);
     for (j=0; j<2; j++) {
       IGABoundary bnd;
       PetscInt    field = 0;

demo/NavierStokesVMS.c

   ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
   ierr = IGAAxisSetOrder(axis0,p[0]);CHKERRQ(ierr);
   ierr = IGAAxisSetPeriodic(axis0,PETSC_TRUE);CHKERRQ(ierr);
-  //ierr = IGAAxisInitUniform(axis0,p,C,N,-0.5*Lx,0.5*Lx);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,p[0],C[0],N[0],0,Lx);CHKERRQ(ierr);
+  //ierr = IGAAxisInitUniform(axis0,N,-0.5*Lx,0.5*Lx);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,N[0],0.0,Lx,C[0]);CHKERRQ(ierr);
   IGAAxis axis1;
   ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
   ierr = IGAAxisSetOrder(axis1,p[1]);CHKERRQ(ierr);
-  //ierr = IGAAxisInitUniform(axis1,p,C,N,-0.5*Ly,0.5*Ly);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis1,p[1],C[1],N[1],0,Ly);CHKERRQ(ierr);
+  //ierr = IGAAxisInitUniform(axis1,N,-0.5*Ly,0.5*Ly);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis1,N[1],0.0,Ly,C[1]);CHKERRQ(ierr);
   IGAAxis axis2;
   ierr = IGAGetAxis(iga,2,&axis2);CHKERRQ(ierr);
   ierr = IGAAxisSetOrder(axis2,p[2]);CHKERRQ(ierr);
   ierr = IGAAxisSetPeriodic(axis2,PETSC_TRUE);CHKERRQ(ierr);
-  //ierr = IGAAxisInitUniform(axis2,p,C,N,-0.5*Lz,0.5*Lz);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis2,p[2],C[2],N[2],0,Lz);CHKERRQ(ierr);
+  //ierr = IGAAxisInitUniform(axis2,N,-0.5*Lz,0.5*Lz);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis2,N[2],0.0,Lz,C[2]);CHKERRQ(ierr);
 
   IGABoundary bnd;
   PetscInt dir=1,side,field;

demo/PatternFormation.c

     ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
     ierr = IGAAxisSetPeriodic(axis,PETSC_TRUE);CHKERRQ(ierr);
     ierr = IGAAxisSetOrder(axis,p[i]);CHKERRQ(ierr);
-    ierr = IGAAxisInitUniform(axis,p[i],C[i],N[i],-1.0,+1.0);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],-1.0,+1.0,C[i]);CHKERRQ(ierr);
   }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   if (n3<3) C[2] = C[0]; if (n3<2) C[1] = C[0];
   for (i=0; i<dim; i++)  if (C[i] ==-1) C[i] = p[i] - 1;
   
-  IGA iga; IGAAxis axis;
+  IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,dim);CHKERRQ(ierr);
   ierr = IGASetDof(iga,dof);CHKERRQ(ierr);
   for (i=0; i<dim; i++) {
+    IGAAxis axis;
     ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
     ierr = IGAAxisSetPeriodic(axis,t[i]);CHKERRQ(ierr);
-    ierr = IGAAxisInitUniform(axis,p[i],C[i],N[i],0.0,1.0);CHKERRQ(ierr);
+    ierr = IGAAxisSetOrder(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],0.0,1.0,C[i]);CHKERRQ(ierr);
     for (j=0; j<2; j++) {
       IGABoundary bnd;
       PetscInt    field = 0;
 
   IGAAxis axis;
   ierr = IGAGetAxis(iga,0,&axis);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis,p,C,N,-1.0,+1.0);CHKERRQ(ierr);
+  ierr = IGAAxisSetOrder(axis,p);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis,N,-1.0,+1.0,C);CHKERRQ(ierr);
 
   IGABoundary left;
   IGABoundary right;
   ierr = IGASetDim(iga,2);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
 
-  IGAAxis axis0;
-  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,p[0],C[0],N[0],-1.0,1.0);CHKERRQ(ierr);
-  IGAAxis axis1;
-  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis1,p[1],C[1],N[1],-1.0,1.0);CHKERRQ(ierr);
-
+  PetscInt i;
+  for (i=0; i<2; i++) {
+    IGAAxis axis;
+    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
+    ierr = IGAAxisSetOrder(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],-1.0,+1.0,C[i]);CHKERRQ(ierr);
+  }
   IGABoundary bnd;
   PetscInt dir,side;
   PetscScalar value = 1.0;
   ierr = IGASetDim(iga,3);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
   
-  IGAAxis axis0;
-  ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,p[0],C[0],N[0],-1.0,1.0);CHKERRQ(ierr);
-  IGAAxis axis1;
-  ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis1,p[1],C[1],N[1],-1.0,1.0);CHKERRQ(ierr);
-  IGAAxis axis2;
-  ierr = IGAGetAxis(iga,2,&axis2);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis2,p[2],C[2],N[2],-1.0,1.0);CHKERRQ(ierr);
+  PetscInt i;
+  for (i=0; i<3; i++) {
+    IGAAxis axis;
+    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
+    ierr = IGAAxisSetOrder(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],-1.0,+1.0,C[i]);CHKERRQ(ierr);
+  }
 
   IGABoundary bnd;
   PetscInt dir,side;
 extern PetscErrorCode IGAAxisSetKnots(IGAAxis axis,PetscInt m,PetscReal U[]);
 extern PetscErrorCode IGAAxisGetKnots(IGAAxis axis,PetscInt *m,PetscReal *U[]);
 extern PetscErrorCode IGAAxisInitBreaks(IGAAxis axis,PetscInt r,PetscReal u[],PetscInt s);
-extern PetscErrorCode IGAAxisInitUniform(IGAAxis axis,PetscInt p,PetscInt C,
-                                         PetscInt E,PetscReal Ui,PetscReal Uf);
+extern PetscErrorCode IGAAxisInitUniform(IGAAxis axis,PetscInt N,PetscReal Ui,PetscReal Uf,PetscInt C);
 extern PetscErrorCode IGAAxisCheck(IGAAxis axis);
 
 struct _n_IGARule {
 
 #undef  __FUNCT__
 #define __FUNCT__ "IGAAxisInitUniform"
-PetscErrorCode IGAAxisInitUniform(IGAAxis axis,PetscInt p,PetscInt C,
-                                  PetscInt E,PetscReal Ui,PetscReal Uf)
+PetscErrorCode IGAAxisInitUniform(IGAAxis axis,PetscInt N,PetscReal Ui,PetscReal Uf,PetscInt C)
 {
   PetscInt       i,j,k;
-  PetscInt       s,n,m;
+  PetscInt       p,s,n,m;
   PetscReal      *U;
   PetscErrorCode ierr;
   PetscFunctionBegin;
   PetscValidPointer(axis,1);
 
-  if (p == PETSC_DEFAULT) p = axis->p;
-  if (C == PETSC_DEFAULT) C = p-1;
-  if (C == PETSC_DECIDE)  C = p-1;
+  if (C == PETSC_DEFAULT) C = axis->p-1;
+  if (C == PETSC_DECIDE)  C = axis->p-1;
 
-  if (E < 1)
+  if (axis->p < 1)
+    SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,
+            "Must call IGAAxisSetOrder() first");
+  if (N < 1)
     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,
-             "Number of elements must be grather than zero, got %D",E);
-  if (p < 1)
-    SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,
-             "Polynomial order must be grather than zero, got %D",p);
-  if (C < 0 || C >= p)
-    SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,
-             "Continuity must be in range [0,%D], got %D",p-1,C);
+             "Number of elements must be grather than zero, got %D",N);
   if (Ui >= Uf)
     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,
              "Initial value %G must be less than final value %G",Ui,Uf);
+  if (C < 0 || C >= axis->p)
+    SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,
+             "Continuity must be in range [0,%D], got %D",axis->p-1,C);
 
+  p = axis->p;
   s = p - C; /* multiplicity */
-  m = 2*(p+1) + (E-1)*s - 1; /* last knot index */
+  m = 2*(p+1) + (N-1)*s - 1; /* last knot index */
   n = m - p - 1; /* last basis function index */
   ierr = PetscMalloc1(m+1,PetscReal,&U);CHKERRQ(ierr);
 
     U[k]   = Ui;
     U[m-k] = Uf;
   }
-  for(i=1; i<=(E-1); i++) { /* (E-1) breaks */
+  for(i=1; i<=(N-1); i++) { /* (N-1) breaks */
     for(j=1; j<=s; j++)     /* s times */
-      U[k++] = Ui + i * ((Uf-Ui)/E);
+      U[k++] = Ui + i * ((Uf-Ui)/N);
   }
   if (axis->periodic) {
     for(k=0; k<=C; k++) { /* periodic part */

test/Test_SNES_2D.c

 
   IGAAxis axis0;
   ierr = IGAGetAxis(iga,0,&axis0);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis0,p,C,N,-1.0,1.0);CHKERRQ(ierr);
+  ierr = IGAAxisSetOrder(axis0,p);CHKERRQ(ierr);
+  ierr = IGAAxisInitUniform(axis0,N,-1.0,1.0,C);CHKERRQ(ierr);
   IGAAxis axis1;
   ierr = IGAGetAxis(iga,1,&axis1);CHKERRQ(ierr);
   ierr = IGAAxisCopy(axis0,axis1);CHKERRQ(ierr);