Commits

Lisandro Dalcin committed 97b0a1c

Update L2 projection 1D/2D demos

Comments (0)

Files changed (2)

   return (x-1)*(x-a)*(x+a)*(x+1);
 }
 
+PetscScalar Hill(PetscReal x)
+{
+  return (x-1)*(x-1)*(x+1)*(x+1);
+}
+
+PetscScalar Sine(PetscReal x)
+{
+  return sin(M_PI*x);
+}
+
 
 typedef struct {
   PetscScalar (*Function)(PetscReal x);
 {
   AppCtx *user = (AppCtx *)ctx;
   PetscReal x = p->point[0];
+  PetscScalar f = user->Function(x);
+
   PetscReal *N = p->shape[0];
-  
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na = N[a];
       PetscReal Nb = N[b];
       K[a*nen+b] = Na * Nb;
     }
-    F[a] = Na * user->Function(x);
+    F[a] = Na * f;
   }
   return 0;
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal x = p->point[0];
+  PetscScalar f = user->Function(x);
+
+  PetscScalar u;
+  IGAPointFormValue(p,U,&u);
+
+  PetscReal e = PetscAbsScalar(u - f);
+  S[0] = e*e;
+
+  return 0;
+}
+
 #undef __FUNCT__
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
   AppCtx user;
 
-  PetscInt choice=2;
-  const char *choicelist[] = {"line", "parabola", "poly3", "poly4", 0};
-  PetscBool t=PETSC_FALSE;
-  PetscInt N=16, p=2, C=PETSC_DECIDE;
-  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Projection1D Options","IGA");CHKERRQ(ierr);
-  ierr = PetscOptionsBool("-periodic", "periodic",__FILE__,t,&t,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-N", "number of elements (along one dimension)",__FILE__,N,&N,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-p", "polynomial order",__FILE__,p,&p,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsInt("-C", "global continuity order",__FILE__,C,&C,PETSC_NULL);CHKERRQ(ierr);
-  ierr = PetscOptionsEList("-function","1D function",__FILE__,choicelist,4,choicelist[choice],&choice,PETSC_NULL);CHKERRQ(ierr);
+  PetscBool w=PETSC_FALSE;
+  PetscInt  N=16;
+  PetscInt  p=2;
+  PetscInt  C=PETSC_DECIDE;
+  PetscInt  choice=2;
+  const char *choicelist[] = {"line", "parabola", "poly3", "poly4", "hill", "sine", 0};
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","L2 Projection 1D Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBool ("-periodic", "periodicity",__FILE__,w,&w,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt  ("-N", "number of elements",__FILE__,N,&N,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt  ("-p", "polynomial order",  __FILE__,p,&p,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsInt  ("-C", "continuity order",  __FILE__,C,&C,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEList("-function","1D function", __FILE__,choicelist,6,choicelist[choice],&choice,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
   if (C == PETSC_DECIDE) C = p-1;
   switch (choice) {
   case 1: user.Function = Parabola; break;
   case 2: user.Function = Poly3;    break;
   case 3: user.Function = Poly4;    break;
+  case 4: user.Function = Hill;     break;
+  case 5: user.Function = Sine;     break;
   }
 
   IGA iga;
   ierr = IGASetDof(iga,1);CHKERRQ(ierr);
   IGAAxis axis;
   ierr = IGAGetAxis(iga,0,&axis);CHKERRQ(ierr);
-  ierr = IGAAxisSetPeriodic(axis,t);CHKERRQ(ierr);
+  ierr = IGAAxisSetPeriodic(axis,w);CHKERRQ(ierr);
   ierr = IGAAxisSetDegree(axis,p);CHKERRQ(ierr);
   ierr = IGAAxisInitUniform(axis,N,-1.0,1.0,C);CHKERRQ(ierr);
 
   ierr = IGACreateVec(iga,&b);CHKERRQ(ierr);
   ierr = IGASetUserSystem(iga,System,&user);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);
-  
+  PetscScalar error = 0;
+  ierr = IGAFormScalar(iga,x,1,&error,Error,&user);CHKERRQ(ierr);
+  error = PetscSqrtReal(PetscRealPart(error));
+  PetscBool print_error = PETSC_FALSE;
+  ierr = PetscOptionsGetBool(0,"-error",&print_error,0);CHKERRQ(ierr);
+  if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);}
+
+  PetscBool draw = PETSC_FALSE;
+  ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
+  if (draw) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);
     /**/ - 1.0/3 * exp(-pow(X+1,2) - pow(Y,2));
 }
 
+PetscScalar Hill(PetscReal x, PetscReal y)
+{
+  return   (x-1)*(x-1)*(x+1)*(x+1)
+    /**/ * (y-1)*(y-1)*(y+1)*(y+1);
+}
+
+PetscScalar Sine(PetscReal x, PetscReal y)
+{
+  return   sin(M_PI*x)
+    /**/ * sin(M_PI*y);
+}
+
 typedef struct {
   PetscScalar (*Function)(PetscReal x, PetscReal y);
 } AppCtx;
   AppCtx *user = (AppCtx *)ctx;
   PetscReal x = p->point[0];
   PetscReal y = p->point[1];
+  PetscScalar f = user->Function(x,y);
+
   PetscReal *N = p->shape[0];
-
   PetscInt a,b,nen=p->nen;
   for (a=0; a<nen; a++) {
     PetscReal Na = N[a];
       PetscReal Nb = N[b];
       K[a*nen+b] = Na * Nb;
     }
-    F[a] = Na * user->Function(x,y);
+    F[a] = Na * f;
   }
   return 0;
 }
 
+#undef  __FUNCT__
+#define __FUNCT__ "Error"
+PetscErrorCode Error(IGAPoint p,const PetscScalar *U,PetscInt n,PetscScalar *S,void *ctx)
+{
+  AppCtx *user = (AppCtx *)ctx;
+  PetscReal x = p->point[0];
+  PetscReal y = p->point[1];
+  PetscScalar f = user->Function(x,y);
+
+  PetscScalar u;
+  IGAPointFormValue(p,U,&u);
+
+  PetscReal e = PetscAbsScalar(u - f);
+  S[0] = e*e;
+
+  return 0;
+}
+
 #undef __FUNCT__
 #define __FUNCT__ "main"
 int main(int argc, char *argv[]) {
 
   AppCtx user;
 
-  PetscInt choice=0;
-  const char *choicelist[] = {"paraboloid", "peaks", 0};
-  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,"","Projection2D 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 = PetscOptionsEList("-function","2D function",__FILE__,choicelist,2,choicelist[choice],&choice,PETSC_NULL);CHKERRQ(ierr);
+  PetscBool flg;
+  PetscBool w[2] = {PETSC_FALSE,PETSC_FALSE}; PetscInt nw = 2;
+  PetscInt  N[2] = {16,16}, nN = 2;
+  PetscInt  p[2] = { 2, 2}, np = 2;
+  PetscInt  C[2] = {-1,-1}, nC = 2;
+  PetscInt  choice=0;
+  const char *choicelist[] = {"paraboloid", "peaks", "hill", "sine", 0};
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","L2 Projection 2D Options","IGA");CHKERRQ(ierr);
+  ierr = PetscOptionsBoolArray("-periodic", "periodicity",    __FILE__,w,&nw,&flg);CHKERRQ(ierr);
+  if (flg && nw==0) w[nw++] = PETSC_TRUE;
+  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","continuity order",  __FILE__,C,&nC,PETSC_NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEList    ("-function","2D function",__FILE__,choicelist,4,choicelist[choice],&choice,PETSC_NULL);CHKERRQ(ierr);
   ierr = PetscOptionsEnd();CHKERRQ(ierr);
+  if (nw == 1) w[1] = w[0];
   if (nN == 1) N[1] = N[0];
   if (np == 1) p[1] = p[0];
   if (nC == 1) C[1] = C[0];
   switch (choice) {
   case 0: user.Function = Paraboloid; break;
   case 1: user.Function = Peaks;      break;
+  case 2: user.Function = Hill;       break;
+  case 3: user.Function = Sine;       break;
   }
 
   IGA iga;
   for (i=0; i<2; i++) {
     IGAAxis axis;
     ierr = IGAGetAxis(iga,i,&axis);CHKERRQ(ierr);
-  ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
-  ierr = IGAAxisInitUniform(axis,N[i],-1.0,1.0,C[i]);CHKERRQ(ierr);
+    ierr = IGAAxisSetPeriodic(axis,w[i]);CHKERRQ(ierr);
+    ierr = IGAAxisSetDegree(axis,p[i]);CHKERRQ(ierr);
+    ierr = IGAAxisInitUniform(axis,N[i],-1.0,1.0,C[i]);CHKERRQ(ierr);
   }
   ierr = IGASetFromOptions(iga);CHKERRQ(ierr);
   ierr = IGASetUp(iga);CHKERRQ(ierr);
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
 
-  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));
+  PetscBool print_error = PETSC_FALSE;
+  ierr = PetscOptionsGetBool(0,"-error",&print_error,0);CHKERRQ(ierr);
+  if (print_error) {ierr = PetscPrintf(PETSC_COMM_WORLD,"L2 error = %G\n",error);CHKERRQ(ierr);}
+
+  PetscBool draw = PETSC_FALSE;
+  ierr = PetscOptionsGetBool(0,"-draw",&draw,0);CHKERRQ(ierr);
+  if (draw) {ierr = VecView(x,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
+
   ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
   ierr = MatDestroy(&A);CHKERRQ(ierr);
   ierr = VecDestroy(&x);CHKERRQ(ierr);