Nathan Collier avatar Nathan Collier committed d469d07

sloppy previous commit

Comments (0)

Files changed (2)

demo/BoundaryIntegral.c

+#include "petiga.h"
+
+#undef  __FUNCT__
+#define __FUNCT__ "System"
+PetscErrorCode System(IGAPoint p,PetscScalar *K,PetscScalar *F,void *ctx)
+{
+  PetscReal *N0 = p->shape[0];
+  PetscReal (*N1)[2] = (PetscReal (*)[2]) p->shape[1];
+  PetscInt a,b,nen=p->nen;
+  for (a=0; a<nen; a++) {
+    PetscReal Na   = N0[a];
+    PetscReal Na_x = N1[a][0];
+    PetscReal Na_y = N1[a][1];
+    for (b=0; b<nen; b++) {
+      PetscReal Nb_x = N1[b][0];
+      PetscReal Nb_y = N1[b][1];
+      K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y;
+    }
+    F[a] = Na * 0.0;
+  }
+  return 0;
+}
+
+#undef  __FUNCT__
+#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;
+  }
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char *argv[]) {
+
+  PetscErrorCode  ierr;
+  ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
+
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson2D Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
+  IGA iga;
+  ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
+  ierr = IGASetDim(iga,2);CHKERRQ(ierr);
+  ierr = IGASetDof(iga,1);CHKERRQ(ierr);
+
+  IGABoundary bnd;
+  ierr = IGAGetBoundary(iga,0,0,&bnd);CHKERRQ(ierr); // u = 0 on [0,:]
+  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
+  ierr = IGAGetBoundary(iga,0,1,&bnd);CHKERRQ(ierr); // grad u . n = h on [1,:]
+  ierr = IGABoundarySetUserSystem(bnd,Neumann,PETSC_NULL);CHKERRQ(ierr);
+  ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
+  ierr = IGASetUp(iga);CHKERRQ(ierr);
+
+  Mat A;
+  Vec x,b;
+  ierr = IGACreateMat(iga,&A);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&x);CHKERRQ(ierr);
+  ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
+  ierr = IGASetUserSystem(iga,System,PETSC_NULL);CHKERRQ(ierr);
+  ierr = IGAComputeSystem(iga,A,b);CHKERRQ(ierr);
+  
+  KSP ksp;
+  ierr = IGACreateKSP(iga,&ksp);CHKERRQ(ierr);
+  ierr = KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
+  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
+  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
+
+  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;
+}
       PetscReal Nb_y = N1[b][1];
       K[a*nen+b] = Na_x*Nb_x + Na_y*Nb_y;
     }
-    F[a] = Na * 0.0;
-  }
-  return 0;
-}
-
-#undef  __FUNCT__
-#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;
   }
   return 0;
   PetscErrorCode  ierr;
   ierr = PetscInitialize(&argc,&argv,0,0);CHKERRQ(ierr);
 
+  PetscInt N[2] = {16,16}, nN = 2; 
+  PetscInt p[2] = { 2, 2}, np = 2;
+  PetscInt C[2] = {-1,-1}, nC = 2;
   ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Poisson2D Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-N","number of elements",     __FILE__,N,&nN,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-p","polynomial order",       __FILE__,p,&np,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsIntArray("-C","global continuity order",__FILE__,C,&nC,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  if (nN == 1) N[1] = N[0];
+  if (np == 1) p[1] = p[0];
+  if (nC == 1) C[1] = C[0];
+  if (C[0] == -1) C[0] = p[0]-1;
+  if (C[1] == -1) C[1] = p[1]-1;
 
   IGA iga;
   ierr = IGACreate(PETSC_COMM_WORLD,&iga);CHKERRQ(ierr);
   ierr = IGASetDim(iga,2);CHKERRQ(ierr);
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
 
+  PetscInt i;
+  for (i=0; i<2; i++) {
+    IGAAxis axis;
+    ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
+    ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],-1.0,+1.0,C[i]);CHKERRQ(ierr);
+  }
   IGABoundary bnd;
-  ierr = IGAGetBoundary(iga,0,0,&bnd);CHKERRQ(ierr); // u = 0 on [0,:]
-  ierr = IGABoundarySetValue(bnd,0,0.0);CHKERRQ(ierr);
-  ierr = IGAGetBoundary(iga,0,1,&bnd);CHKERRQ(ierr); // grad u . n = h on [1,:]
-  ierr = IGABoundarySetUserSystem(bnd,Neumann,PETSC_NULL);CHKERRQ(ierr);
+  PetscInt dir,side;
+  PetscScalar value = 1.0;
+  for (dir=0; dir<2; dir++) {
+    for (side=0; side<2; side++) {
+      ierr = IGAGetBoundary(iga,dir,side,&bnd);CHKERRQ(ierr);
+      ierr = IGABoundarySetValue(bnd,0,value);CHKERRQ(ierr);
+    }
+  }
+
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);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.